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Introduction 


This is a request for a second renewal (3 rd year of funding) of a research project on the 
topic of multiresolution and explicit methods for vector field analysis and visualization. 
In this report, we describe the progress made on this research project during the second 
year and give a statement of the planned research for the third year. 

Report on Research Progress During the Second Year 

There are two aspects to this research project. The first is concerned with the 
development of techniques for computing tangent curves for use in visualizing flow 
fields. The second aspect of the research project is concerned with the development of 
multiresolution methods for curvilinear grids and their use as tools for visualization, 
analysis and archiving of flow data. We report on our work on the development of 
numerical methods for tangent curve computation first. 

The basic idea of these methods is to decompose the domain into a collection of 
triangles (or tetrahedra) and assume linear variation of the vector field over each cell. 
With this assumption, the equations which define a tangent curve become a system of 
linear, constant coefficient ODE’s which can be solved explicitly. In our report on the 
progress made in the first year of the project, we reported on our work in 2D. We have 
done the mathematical analysis on the explicit and incremental methods for 2D and 
written and tested thoroughly the computer programs implementing certain aspects of 
these algorithms. We have submitted for publication (see [1]) some of our research 
results in this area. We mentioned in the previous progress report our plans to extend our 
work in this area to 3D. This involves tetrahedral decompositions of the 3D curvilinear 
grids (see our submitted research in this area [2] ) and the development of explicit and 
incremental methods for linearly varying vector fields over tetrahedra. We have 
completed the mathematical analysis and are in the program development stages at this 
moment. We include here the mathematics for the 3D case. In the situation of a 2D flow 
over a 2D domain, there were five (5) distinct cases determined by the eigenvalues of the 
coeficient matrix of the ODE defining the tangent curves of the flow. In the 3D situation 
there are nine (9) cases. 
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Case 4) A has one real, nonzero and two complex, eigenvalues (r, * 0, p ± Xi, X * o) 
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Case 7) A has two real, nonzero and one zero eigenvalues (r, = 0, 0 * r 2 ^ r 3 * o) 


P(t) = E x t + E 2 (e' r > - 1) + £ 3 (e' r ' -l) + P 0 



Case 8) A has two real equal, nonzero and one zero eigenvalues (r, =0, r 2 = r z ^ o) 

P(t) = E x t + E 2 (e ,r ' - 1) + E i te"’ 1 + P 0 




Based upon the explicit representations give above, we have developed incremental 
methods which compute points on the exact solution of the tangent curves. We present 
the formulas for these incremental methods for the nine (9) cases below. (The matrix, E, 
consists of columns which are the eigenvectors of matrix A.) 
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3D Incremental Method in Cartesian Coordinates 
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Case 5) A has one real, nonzero and two zero eigenvalues (r, = r 2 = 0, r 3 * O) 
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Case 9) A has one zero and two complex, eigenvalues (r, = 0, p ± Xi, X * 0) 
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We did not discuss this in the report on the first year's progress, but we have also 
developed incremental methods that compute points on the true tangent curves (not 
approximations) using barycentric coordinates. The reason for doing the computations in 
barycentric coordinates is to make it cheap and easy to determine when the curve leaves 
one triangle and enters a neighboring triangle. All that needs to be done is to check the 
sign of the barycentric coordinates. If all three are not positive, then the curve has left the 
triangle and which barycentric is not positive tells us which neighboring triangle the 
tangent curve has entered. We give the formulas for the 2D incremental methods using 
barycentric coordinates below. These formulas are based upon explicit, exact solutions 
represented in terms of barycentric coordinates which we also include here. 

In the formulas below, p(t) represents the tangent curve in barycentric coordinates. 
That is 


p( 0 


Uiit)' 
bj{ t) 
Vb k {t)J 


and p{t) = bi(t)Vj + bj(t)Vj + b k (t)V k 


where V i ,V j ,V k are the velocity vectors at the vertices of the triangle currently 

containing the tangent curve. The matrix V (used below) includes the information which 
defines the tangent curve which is similar in role to the matrix A and column vector B in 
the cartesian coordinate case. 


2D Explicit Solutions in Barycentric Coordinates: 

p'(t)=Vp(t), p(0) = p 0 

Case 1) V has two real, nonzero and one zero eigenvalues 0 * r, * r 2 * 0, r 3 = 0 
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We are at this moment in the process of finishing off the development of the 
equations similar to the above for 3D and plan on implementing all of the 3D methods 
mentioned above in the near future. 
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We now move to a discussion on our research progress in the other topic of this 
research project and this is in the area of multiresolution modeling. Our work in this area 
can be broken down into two subareas: i) triangular grids with applications to arbitrary 
domains and ii) curvilinear grids. We first take up the topic of triangular grids and our 
current extensions to tetrahedral grids. 

Earlier we reported on some new Haar type wavelets for subdivision of triangular 
domains as illustrated in Figure 1. The basic decomposition step is illustrated in Figure 2. 
The fs are the characteristic functions for the subtriangles and therefore form a basis for 
piecewise contant functions. The vp's are the wavelet functions which capture the detail 
(error) and hopefully they are mutually orthogonal so that we obtain best least squares 
approximations when we do the wavelet decompositions. 





Figure 1 . Nested triangular domain with labeling scheme. 
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Figure 2. Basic decomposition step for Haar wavelets over nested triangular domains 

Biorthogonal wavelets within this context have previously appeared. We have 
developed some improvements. While these new wavelets are also only biorthogonal, 
they have the property that when the areas do become the same, the wavelets become 
mutually orthogonal and this is a definite plus when it comes to ability to approximate 
data. A nice example of where this happens is with the standard triangular decomposition 
of the sphere. 

We have applied our new wavelets to data consisting of a vector field defined over the 
sphere. Our research on Haar wavelets for nested triangular grids and applications to 
fluid flow analysis over a spherical domain is reported in the submitted manuscript [4] . 



View 1,1% View 1, 20% 

Figure 3. Partial (nearly orthogonal) spherical wavelet reconstruction of a flow field 

defined over a spherical domain. 
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Also we have applied our new, nearly orthogonal wavelets to some real data 
provided to us by Roger Crawfis and Nelson Max of LLNL. This is simulated data 
resulting from a global weather model. A typical image for this work is shown in Figure 
4. More details can be found in the submitted manuscript [4], 



View 1,100% View 1,3% 

Figure 4. Partial reconstruction of global weather data. 


We have started our work on extending the research on wavelets defined over nested 
triangular grids to the case of nested tetrahedral grids. One surprise has already come up. 
We had thought that it would be natural to use the nested tetrahedral grids based upon the 
basic decompositon as shown in Figure 5. This decompostion seemed to be the 3D 
analog of the decompositon used for triangles shown in Figure 1. We started to work 
with these types of grids and then found out much to our surprise that we could not define 
affine invariant wavelets over these types of decompositions. We feel that the property of 
affine invariance is extremely important. In essence it says that the results of any analysis 
or compression based upon these wavelets is not affected by the way we label the 
triangles. Clearly, it is undesirable to have an algorithm that depends on the inner details 
of how it is coded. In order to maintain the property of affine invariance, we have moved 
to a different type of basic decomposition step. It is shown in Figure 6. Rather than 8 
subtetraheda per tetrahedron, we have 17 (yes 17!) subtetrahedra per tetrahedron. 
Utilizing Mathemaica we have recently been able to define some affine invariant Haar 
wavelets based upon this type of decomposition. We are really quite exicited about these 
results and can’t wait to implement them and see how they perform. 
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Figure 6. Decomposing a tetrahedron into 17 sub-tetrahedra. The top figure shows the 
edges on the faces and the bottom is an exploded view. 

In order to convey some idea of our current progress on Haar wavelets defined over 
tetrahedrizations shown in Figure 6, we include some of our preliminary results for the 
basic refinement equations for these new wavelets. The free parameters below are c, ao, 
ai, a 2 , b 0 , bi, b 2 , and b 3 . We impose the orthogonality conditions and then use 
Mathematica to solve for the refinement equations with the values give below. 
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Solutions: 

The following constants are used in the solutions which give values for ao, &i> a 2 , bo, b], 


b 2 , and b 3 . 


Constant 

Symbolic Value 

ki 

— ±— Vl7 
16 64 

k2 

— ±— Vl7 ±-^a/221 
16 12 48 

k 3 

_C ± 13C Vl7 
16 192 

k4 

— ±-Vl7 
16 16 


Solution 1: 

-7/1 Oc - 51/5 k! 
ki 

-21/80c - 16/5 k! 
-3/40c - 1/5 k! 
-3/40c - 1/5 k! 
211/80c + 216/5 ki 
-13/80c - 8/5 ki 
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bo 

bi 

b2 

b 3 
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Solution 2: 
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Solution 3: 
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b 3 
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136c 2 


We now discuss our research in the area of multiresolution models for curvilinear 
grids. In the original proposal we mentioned the idea of using multiresolution models for 
parametric curves in this context. The idea is based upon the observation that a 
curvilinear grid can be viewed as a parametric map and researchers have already 
developed techniques for using wavelets on parametric curves. Unfortunately there is a 
basic problem with this approach and that lies in the fact that not only do you obtain a 
lower resolution approximation to the flow, but you also get a lower resolution 


17 




approximation to the boundary of the domain. This means that the domam will be 
"smoothed down" also. We have overcome this problem in 2D by using a knot removal 
technique to yield a collection of nested domains as illustrated in Figure 7. (Note that the 
detail of the inner boundary (not the outer, though) remains through all levels of 
approximation) . 





Figure 8. Decomposition by knot removal 

We have developed some Haar wavelets for these types of decomposition and 
examples of our new algorithms are included in the submitted manuscripts [1] and [4], 
We plan to write another research paper delineating more of the details of these methods, 
but will most likely wait to do this after our 3D work is completed. We have done most 
of the developmental work for the 3D case and are in the middle of doing the 
implementation. 

As we had planned, we have extended the Haar wavelets to higher degree continuity. 
We have used the so called "lifting scheme." These higher order methods are more 
expensive to compute, but they are much more efficient. An example is shown in Figure 

9. 
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Figure 9b. Improved linear wavelet approximations over nested curvilinear grid. 
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Summary of 2 nd Year Research 


• Explicit Methods 

We have developed, implemented and thoroughly tested the 2D explicit 
techniques based upon linear variation over triangles. We can compute a single tangent 
curve or compute a complete topological graph of a 2D flow over a curvilinear grid. 

We have done the mathematical analysis for the extension of the explicit method 
to 3D (some of the equations are given earlier in this report) and we are currently 
developing and implementing incremental methods for 3D tetrahedra and curvilinear 
grids. 

• Multiresolution Methods for Curvilinear Grids 

We extended our piecewise constant wavelets over 2D curvilinear grids to 
piecewise linear wavelets. We are currently in the middle of implementing and testing 
these methods. Some preliminary results are shown in Figure 9 of this report. 

We have developed the data structures and algorithms for extending our knot 
removal algorithms for curvilinear grids to 3D. This requires the tetrahedrization of 
general polyhedra with fixed boundary polygons. We are now in the middle of 
implementing these results. 

• Multiresolution methods for nested triangular grids 

We have used a lifting scheme to extend the 2D Haar wavelets to piecewise 

linear. 


We have developed the basic decompositionn and refinement equations for Haar 
3D wavelets over certain (see Figure 6) tetrahedral subdivisions. These new equations 
appear for the first time in this report. 
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Publications resulting from the research of this project 

(All publications listed here specifically acknowledge support under this NASA grant and 

copies are attached) 

1. G. M. Nielson, I.-H. Jung, N. Srinivasan, J. Sung and J.-B. Yoon, Tools for 
computing tangent curves and topological graphs for visualizing piecewise linearly 
varying vector fields over triangulated domains, to appear in the book, Scientific 
Visualization: Overviews, Methodologies, and Techniques, CS Press, 1997. pp. 517- 
558. 

2. G. M. Nielson, Tools for triangulations and tetrahedrizations and construction 
functions defined over them, to appear in the book. Scientific Visualization: 
Overviews, Methodologies, and Techniques, CS Press, 1997, pp. 419-515. 

3. G. M. Nielson and Junwon Sung, Interval volume tetrahedrization, submitted to 
Visualization '97 Conference. 

4. G. M. Nielson, Il-Hong Jung and Junwon Sung, Haar wavelets over triangular 
domains with applications to multiresolution models for flow over a sphere, 
submitted to Visualization '97 Conference. 

5. G. M. Nielson and Richard Franke, Computing the separating surface for segmented 
data, submitted to Visualization '97 Conference. 
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Statement of Work for Next Year 

As we have previously done, we structure our discussion upon the topics of this 
research proposal: explicit methods, multiresolution methods for curvilinear grids, and 
multiresolution methods for triangular grids. 

As far as explicit methods go, in the past year, we have created some algorithms for 
the efficient computation of 2D tangent curves based upon our ideas of explicit solutions 
to linearly varying vector fields. One of the types of methods we have developed are 
(exact) incremental methods. We have done all of the mathematical work for the 
extension of these methods to 3D and plan to implement and test these new 3D methods. 
In addition, we have been testing some 2D methods which work directly in barycentric 
coordinates. Barycentric coordinates allow one to discuss tangent curves (and other 
descriptive visualization objects) in a context that is independent of the coordinate system 
used for the data. This has several advantages including the ease of detecting which cell a 
point (on a tangent curve) lines within. We are at this moment in the middle of doing the 
mathematical development for the 3D situation. This is taking some effort, but it appears 
that everything will eventually work out. Once the equations have been obtained, we will 
implement the methods and compare and test them. 

Points of attachment and detachment are special critical points (zero flow) on the 
interior boundary (car body or airplane wing for example). They are very important in the 
computation of topological methods. In the past, most of the literature on this topic was 
not so rigorous. Using the exact solutions for linearly varying vector fields, we have been 
able to give very precise definitions which lead to effective algorithms in the case of 2D 
flows. We reported on this work in [1]. We are planning to extend this work to 3D. 
Here, points of attachment and detachment will be replaced with edges of attachment and 
detachment. We are hopeful that the same rigorous characterization and resulting 
effective algorithms well exist in the 3D case. 

Earlier, we developed some knot removal techniques which lead to a collection of 
nested domains for 2D curvilinear grids. These nested domains were the necessary for 
the multiresolution (wavelet) analysis for vector fields defined over 2D curvilinear grids. 
We subsequently developed, implemented and reported on this research. See [1] and [4]. 
We want to extend all of this to 3D. We have already done the development for the 3D 
nested grids and we now currently extending the Haar wavelets to this context. There 
have been several snags along the way, but we feel that we can overcome this obstacles 
and obtain the results which are comparable to our 2D work. We will implement these 
new 3D Haar wavelets and test them on real data. Just as we were able to use "lifting" 
techniques to extend the Haar wavelets to linear wavelets for the nested 2D curvilinear 
grids, we hope to do the same for 3D grids. 

Previously, we developed some new Haar wavelets for nested triangular domains as 
shown in Figure 1. We have applied these new techniques to a variety of data sets 
including slices of 3D flow around and automobile and weather data representing flow 
over the earth. We have reported on these results in the submitted manuscripts [1] and 
[4], We are currently involved in testing the extension of these methods to higher order, 
linear methods. For the next year, we want to extend all of these results to 3D. We were 
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originally frustrated with the basic mathematics when we first attempted to use the 
decomposition of Figure 5. This did not work out, but we are currently having success 
with the decomposition shown in Figure 6. Some of the early mathematical results are 
mentioned above. We will also use the basic building blocks of these techniques to 
develop some methods for the adaptive volume modeling of 3D volume data. 
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Summary of Proposal for Research in 3 rd Year 


• Explicit Methods 

Implement and test algorithms for tetrahedra based incremental methods. 

Continue development and implementation of 3D methods utilizing Barycentric 

coordinates. 

Develop the proper definitions and characterizations of the "edges of attachment" 
and "edges of detachment" for 3D flows near an inner boundary. 

• Multiresolution Methods for Curvilinear Grids 

We have previously developed the 3D knot removal techniques for generating 
nested domains for 3D curvilinear grids. We are also now doing the development of the 
3D Haar wavelet in this context. We plan to complete this development and implement 
and test these algorithms. 

We propose to extend Haar wavelets for 3D curvilinear grids to piecewise linear 
wavelets for 3D curvilinear grids. 

• Multiresolution methods for nested tetrahedral grids 

Using the subdivision as shown in Figure, we plan to further develop, implement 
and test Haar wavelets for tetrahedral grids. We plan to extend these techniques for 
higher order (linear) wavelets which should have more efficient compression and 
approximation capabilities. 

We plan to use the basic multiresolution techniques to develop some adaptive 
methods for modeling volume data. This would allow details to be added and removed as 
the user's attention is focused in different regions of the data. 
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Budget 

We are requesting funds for the principal investigator as follows: 

Gregory M. Nielson 

2 summer months 

(see budget sheet for amounts) 

We are requesting funds for two students: 

Two ASU Graduate Students 

9 months half time, 3 months full time 
(see budget sheet for amounts) 

We are requesting funds in the operations category as explained on the budget sheet. 
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ESTIMATED BUDGET 


PERIOD OF PERFORMANCE: January 1,1995 through December 31, 1997 








Year 1 

Year 2 

Year 3 

Summary 

I. 

SALARIES & WAGES 










A. PI: Nielson 

0.00 

pm/AY + 

2.00 

pm/Sum 

$16,610 

$17,607 

SI 8,663 

552,880 


B. Co-PI: 

0.00 

pm/AY + 

0.00 

pm/Sum 

$0 

SO 

so 

so 


C. Co-PI: 

0.00 

pm/AY + 

0.00 

pm/Sum 

$0 

$0 

$0 

so 


D. Co-PI: 

0.00 

pm/AY + 

0.00 

pm/Sum 

$0 

$0 

so 

so 


E. Postdocs 0 @ 

0.00 

pm/CY 



SO 

SO 

SO 

so 


F. Technicians 0 @ 

0.00 

pm/CY 



$0 

$0 

SO 

so 


G. Graduate Associates 2 @ 

9.00 

pm/AY + 

6.00 

pm/Sum 

535,020 

$36,050 

$37,767 

5108,837 


H. Graduate Assistants 0@ 

0.00 

pm/AY + 

0.00 

pm/Sum 

$0 

$0 

SO 

SO 


I. Undergraduates 0 @ 

0.00 

pm/AY + 

0.00 

pm/Sum 

SO 

$0 

SO 

so 


J. Secretary 0 @ 

0.00 

pm/CY 



$0 

so 

so 

SO 


TOTAL 





$51,630 

$53,657 

556,430 

S 1 6 1 ,717 

II. 

FRINGE BENEFITS 










A. Faculty at 

25% 




$4,153 

$4,402 

$4,666 

$13,221 


B. Staff at 

30% 




SO 

$0 

so 

SO 


C. Students at 

3% 




$1,051 

SI, 082 

$1,133 

S3, 266 


TOTAL 





$5,204 

$5,484 

$5,799 

S16.487 

III. 

TRAVEL 





$0 

SO 

SO 

SO 

IV. 

OPERATIONS 










A. Consultants 





$0 

$0 

SO 

so 


B. Publication/Page Charges 





$0 

$0 

SO 

SO 


C. Photocopy, telephone, postage, 

mi sc. 




S500 

$500 

S500 

51,500 


D. Material & Lab Supplies/Equipment Use Fees 



SI, 000 

SI, 000 

S 1 ,000 

53,000 


E. Computer Software 





S200 

SO 

SO 

S200 


F. Lab Equipment under $500 





$0 

SO 

SO 

SO 


G. Subcontracts 





$0 

$0 

SO 

so 


H. Participant Support Costs 





$0 

$0 

so 

so 


I. Tuition 0 @ 

$1,895 

per AY 

(4% escalation) 

SO 

so 

so 

so 


J. Student Stipends 





$0 

so 

so 

so 


K. Other 





$0 

so 

so 

so 


TOTAL 





$1,700 

SI, 500 

$1,500 

S4.700 

V. 

CAPITAL EQUIPMENT 





so 

so 

so 

SO 

VI. 

TOTAL DIRECT COSTS 





$58,534 

$60,641 

563,729 

5182,904 

VII. 

INDIRECT COSTS 








S95 ,733 




Year l 

52.0% 

MTDC 

S30.438 







Year 2 

52.5% 

MTDC 


S31,837 






Year 3 

52.5% 

MTDC 



$33,458 


VIII. 

TOTAL PROJECT COSTS 





$88,972 

$92,478 

S97.187 

$278,637 




Chapter 21 

Tools for Computing Tangent 
Curves and Topological 
Graphs for Visualizing 
Piecewise Linearly Varying 
Vector Fields over 
Triangulated Domains 


Gregory M. Nielson, ll-Hong Jung, Nat Srinivasan, Junwon Sung, 
and Jong-Beum Yoon 


Abstract. We describe some methods for computing tangent curves and topological 
graphs for a vector field defined over a two-dimensional domain. We assume that the 
vector field is piecewise linearly defined over a triangulation of the domain. Piecewise 
explicit representations in terms of elementary transcendental functions form the basis of 
algorithms for determining and displaying tangent curves. Topological methods which link 
critical values with separating tangent curves are developed and discussed. 


21.1 Introduction 

Streamlines are a well-established visualization technique for investigating a vector field. 
In this chapter we describe some new techniques for computing these invariant tangent 
curves. Conventional methods for computing tangent curves consist of using numerical 
methods for solving vector- valued initial valued problems. Euler s method is the simplest 
and Runge-Kutta-type methods are often used in practice (see [3,24,39,41,49]). The issues 
involved in selecting a particular algorithm usually focus on the two, often opposing, re- 
quirements of accuracy and speed. Accuracy is especially important in the computation of 
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tangent curves for topological methods. Erroneous results can easily occur unless special 
provisions are taken to control errors. Speed is particularly important for interactive meth- 
ods, but accuracy should not be discounted too much in this context since the goal of the 
visualization is to impart meaningful and valid information. The methods described in this 
chapter allow for a reasonable balance to be achieved between accuracy and speed. 

A streamline (or tangent) curve is a parametric curve 




that is everywhere tangent to the vector field. If 

V{x,y) = 

represents a static vector field, then P is characterized by the equations 


u(x, y) 
v(x,y) 


P'(t) = V(P(t)) = ( 


) ( v[x (t),y(t)) 


(*( 0 , 2 /(<)) 


( 21 . 1 ) 


"Typically there is an entire family of solutions for Equation (21.1) and a particular solution 
is selected with the initial condition 


P(0) 


■c(O) 

2 /( 0 ) 


= Po 


Xq 

yo 


In this chapter, we discuss methods for the special case where V(x, y) is a piecewise 
linear function. The domain is decomposed into a collection of triangles and V has the 
form 


V(x,y) 


u(*,t/) \ = / aiix + awy + bi \ _ V x \ 
v{x,y) ) V a 2 \x + a 2 2 y + 62 )~ \ y ) 


( 21 . 2 ) 


over each triangle. With this assumption, the tangent curve becomes a piecewise concate- 
nation of curve segments with each segment associated with a particular triangle. The entry 
point on a particular triangle provides the initial conditions for the constant coefficient ODE 
which characterizes the curve segment for each triangle. The exit point serves as the entry 
point for the next adjoining triangle domain. This basic idea is further illustrated in Fig- 
ure 21.1. In this way, it is possible to completely characterize and know a tangent curve by 
a sequence of entry (exit) points. 

There are a variety of sources for the type of data covered here. If a 2D flow (or a 3D 
flow assumed to be constant in one direction) is measured at a collection of scattered, planar 
points, then the domain points can serve as the vertices for a triangulation of the domain. 
The Delaunay triangulation of the convex hull would be a possibility. See Nielson [33]. We 
will show examples later of a flow over a spherical domain which represents wind speed 
and direction over the Earth. If the data is measured at scattered locations, these points can 
serve as the vertices of a triangulation of the sphere and the methods of this chapter can be 
applied. If the data is modeled, then the model can be evaluated on a triangular grid and 
again the methods of this chapter will apply. Curvilinear grids which normally associate 
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Figure 21.1: A composite tangent curve running across three triangles. 


with simulated data can also be triangulated by simply inserting one diagonal or the other 
into the quadrilateral cells. See Figure 21.2. 

The values A and B of Equation (21 .2) for a particular triangle are determined from 
the values of the vector field at the vertices of the triangle. If the vertices of the triangle are 
labeled as in Figure 2 1 .3, then the equations 


an x i + a 12Vi + &1 

= u{x it yi) 

an xj + ai22/j + 

£ 

II 

anXk H- ai2 yk + 

= u(x k ,y k ) 

anXi + a 2 2Vi + b 2 

£ 

II 

a 2 i Xj + a 2 2 Vj + ^2 

= v{xj,yj) 

&2l x k + <^22Vk + ^2 

= v(x k ,y k ) 


will yield the values of A and B. It is also possible to use barycentric coordinates to find 
these values. The details of this will be discussed later in Section 2 1 .4. 

The remainder of the chapter is organized as follows. In Section 21.2, we discuss 
explicit methods which are based upon the fact that the differential equation which char- 
acterizes a tangent curve is first order with constant coefficients and so an explicit solution 
in terms of common transcendental functions is possible. Section 2 1 .3 is concerned with 
incremental methods which produce a sequence of points on the curve (or approximations) 
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Figure 21.2: Triangulated curvilinear grid. 


with the computation of each subsequent point based upon the previous point. Section 2 1 .4 
covers much of the same material of Sections 21.2 and 21.3, but within the context of 
barycentric coordinates rather than the conventional Cartesian coordinates. Topological 
methods are covered in Section 21.5. 


21.2 Explicit Methods 


In this section we describe what we call ‘explicit” methods. Because of the special nature 
of the ODE which characterizes a tangent curve, it is possible to obtain explicitly defined 
solutions and therefore, there is no need to use numerical methods to compute the values of 
these curves. Library routines for computing common transcendental functions can be used 
instead. Computing the exit point of a particular curve segment, which will serve as the 
initial point for the next curve segment, amounts to computing the intersection of a curve 
and a line segment. This is equivalent to a univariate root computation problem. We discuss 
two general approaches to this intersection point calculation problem. Each approach is 
covered in a subsequent subsection. In Section 21.2.1 we discuss the approach which is 
based upon a parametric representation of the tangent curve and an implicit representation 
of a line containing one of the edges of the triangle. The parametric curve is substituted 
into the implicit line equation yielding a single equation in the parameter of the tangent 
curve. The root of this equation yields the parameter value of the intersection point and 
this point is subsequently tested to determine if it is actually on the edge (a subset of the 
line). In Section 21.2.2, the second general approach is covered. It is based upon an 
implicit representation of the tangent curve and a parametric representation of the edge of 
the triangle. In this approach, the parametric representation of the edge is substituted into 
the implicit representation of the tangent curve, yielding again a single, univariate equation. 
The root is computed and used to evaluate the parametric representation of the edge so as 
to obtain the intersection point. 
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( u (Wj)> ^Wj)) 


Figure 21.3: Notation and conventions used for data. 


21.2.1 Parametric Tangent Curve, Implicit Edge 


In this section, we discuss the general properties of a tangent curve for a linearly varying 
vector field over a single triangular domain. In particular we give the details of a parametric 
representation, 

«*>-($)• 


(Implicit representations are covered in the next section. Section 21.2.2.) The curve enters 
the triangle at a point P(0) = P 0 and then either attaches to a critical point in the triangle or 
exits from the triangle at a point P(t e ) = P e . We assume that the entry point is given and 
so we need to compute, t e , P e , and the various coefficients of a parametric representation 
of the curve from P 0 to P e . The value t e along with the parametric representation are used 
for displaying the curve and the value P e serves as the entry point for the next triangle. If 
P e lies on one of the edges, then t e must satisfy one of the equations 


fa(x{te),y{te)) = 0. a = i,j,k (21.4) 


where 


fi(x,y) = (x - Xj)(y k - Vj) ~ {y ~ yj)( x k ~ ^j) 

fj(x,y) = (x -**)(»< -yk) ~ {y~ Vk){xi -Xk) 
fk(x,y) = {x - Xi)(yj - yi) - (y - yi)( x j ~ x i)- 

But not all of the roots of f a (x(t ) , y(t )), a = i, j, k are candidates to be t e . It must be the 
case that t e is the smallest of these roots and that P e is actually on the appropriate edge 
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and not simply on the line containing this edge. The point P e could potentially be on any 
of the edges, including the edge containing P 0 , but some simple tests and modifications 
can possibly eliminate some edges and portions of edges where P e might be found. Along 
an edge, the flow is either always in, always out, or there is a special point where the 
flow is parallel to the edge and the direction of flow relative to the triangle changes at this 
point. This is further illustrated in Figure 21.4. The change of direction point P& is easy to 
compute. For example, if the direction changes on edge e k , then 

Pa = tPi + (1 - t)P : (21.5) 

where 


tVi + (1 - t)Vj = e(Pj - P t ) 
for some constant, c, and 0 < t < 1. 




P. 




No edges to search Only a subset to search 


Figure 21.4: Exit point search strategy regions. 
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We now turn our attention to the details of the parametric representation of the tangent 
curve. The tangent curve 



over a particular triangular cell is defined as 



( 21 . 6 ) 


with the initial conditions 



The 2x2 matrix A and the 2 x 1 vector B are determined as in Equation (2 1 .3) or equivalent 
other means. (For example, see Section 21.4 on barycentric coordinates.) The general 
solution of this initial value problem is of the form 


( y(t) ) = P ® = ^ l( ' t)El + + c (2L7) 

where the particular functions <$i, and the coefficients E\, E 2 , and C depend on the 
eigenvalues of A. There are five separate cases: 


Case 1. A has two real, nonzero eigenvalues, 0 ^ r x ^ r 2 7 ^ 0. 

P(t) = e tri E\ + e tr *Ei + P Cl (21.8) 


* = &%)<*- FA Bi = (^rr) (P “ " p ‘ h AP < +B =° 

Case 2. A has one zero and one nonzero eigenvalue, 0 = r\, r 2 ^ 0 . 

P(t) = tEi + ( e tr 3 - 1) E 2 + Po, ( 21 - 9 ) 



Case 3. A has only a zero eigenvalue, r 1 = 0 = r 2 . 

P(t) =tEi + ( 21 - l °) 


Ei = APq + By E-} — g * 
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Case 4. A has a single real, nonzero eigenvalue, rq = r 2 / 0. 

P(t) =e tr 'E l +te tr 'E 2 + P c , (21.11) 

E 1 = P 0 -P C , E 2 = (A - r x I) (P 0 - P c ) , AP c + B = 0. 

Case 5. A has complex eigenvalues, n + Xi, n - X i, A ^ 0. 

P[t) — cos(Xt)e ,it Ei + sin(Xt)e >,t E 2 + P c (21.12) 


Ei = Pa -Pc , E, = (Po - Pc) , AP C + B = 0. 

In Cases 1, 4, and 5, .4 is nonsingular and it is guaranteed that there is a unique critical 
value satisfying AP - + B = 0. In Cases 2 and 3, it is possible that no critical values exist 
or even that an entire line of critical values exist. Not all cases occur with equal frequency. 
Cases 1 and 5 are the predominate ones, Cases 2 and 4 are much less likely to occur, and 
Case 3 is extremely rare. This is explained by taking a look at Figure 2 1 .5. In this example, 
we determined the case for a triangle where we varied the value of the vector field at one 
vertex of the triangle. The flow at two of the vertices is fixed and illustrated by the arrows 
drawn at these vertices. The flow at the vertex marked with the white circle is taken be a 
vector which is based at the vertex and emanating out to an arbitrary point in the plane. 
We classify this point on the basis of the case classification associated with this flow. The 
interior to the parabolic bounded region is Case 5. The boundary is Case 4 except for the 
degenerate subcase of Case 3 indicated by the black box on the boundary. Case 2 consists 
of a line tangent at this point. Of course the tolerances used for determining these regions 
affect their relative occurrence in actual practice. It is interesting and instructive to look 
at some examples of tangent curves for the various cases. We have included samples for 
each of the cases in Figure 2 1 .6. The data for these examples is the same as in Figure 2 1 .5. 
The boxes of Figure 21.5 indicate the tip of the one vector of the triangle for the particular 
example of Figure 21.6. In each of these images, the triangular domain is shown along 
with arrows at each vertex indicating the flow at these points. Particular tangent curves are 
determined by using initial values along a particular edge. We use 10 equally spaced points 
along these edges. The flow is shown at each of these initial points and the tangent curve 
is traversed out in both negative and positive parameter domain for a fixed amount. Some 
additional examples are shown in Figure 21.7. For these examples, the flow data at each 
vertex is indicated by the line segment and again we use initial values at 10 equally spaced 
values along an edge and the tangent curves are traversed out in both positive and negative 
parameter directions. 

Displaying the tangent curve and computing the roots for the determination of t e de- 
pends upon the evaluation of P[t) given in its various forms by the five cases above. De- 
pending upon the computing resources available, it is possible to structure efficient meth- 
ods for performing these computations. In some cases, it is desirable to compute P(iAt), 
i = 0,1,2,..., where At is some fixed increment of the parameter t. This may be the case 
in a situation where these values are used to display the curve and the curve is determined 
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to exit the triangle once one of the functions / a , /j, or fk (of Equation (21.4)) is found 
to change sign from aAt to ( a + 1)A t. The exponentials of Cases 1, 2, 4, and 5 can be 
computed by the computation of a single exp(A<) and subsequent multiplications. The sin 
and cos functions of Case 5 can be computed with the formulas 

sin(t + At) = cos(t)sin(At) + sin(t)cos(At) 

cos[t + At) = cos(t)cos(At) - sin(t)sin(At). 

We should also point out that the formulas which form the basis of the incremental methods 
of Section 2 1 .3.2 can also be used to evaluate the curve at equally spaced parameter values 
for display purposes. 

21.2.2 Parametric Edge, Implicit Tangent Curve 

As we have previously mentioned, there are two general approaches to computing the exit 
point P(t e )- A parametric representation for the tangent curve can be substituted into 
an implicit representation for the line segment as in Section 21.2.1. Or a parametric (or 
explicit) representation of the line segment can be substituted into an implicit representation 
of the curve. In this section, we discuss the latter approach. In this approach, either an 
explicit, y = ax + b or x = cy + d, or a parametric (x(s), y(s)) = sPi + (1 - s)Pj , is used 
to represent the line segment and an implicit representation 

F(x(t),y{t)) = 0 
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Figure 21.7: Additional examples of tangent curves. (Upper-left) Case 1 with n • r 2 < 0. 
(Upper-right) Case 1 with r\ ♦ 7*2 > 0. (Lower) Case 2. 
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is required for the curve. For the edge e k , this leads to one of the following root finding 
problems: 6 

i) F(x,ax + b), min {xi, xj} < x < max{ x x , } 

ii) F(cy + d,y), rnin{y ityj } < y < max{y uyj } 

iii) F(x(s),y(s)), 0 < s < 1. 

Up to now, the missing component for this approach is an implicit representation of the 
streamline curve. We now take up this topic. Again, there are five separate cases. In each of 
the five cases covered below, there are specified two vectors, and E 2 . We let E = (Ei, 

E 2 ). If \E\ ^ 0 then the two vectors Ei and E 2 are linearly independent, and we can use 
the change of variables 



The point C is either P 0 for Cases 2 and 3 or P c for Cases 1, 4, and 5. We now give the 
details for each separate case. 

Case 1. (0 ^ r x ^ r 2 ^ 0) 

If | £~ | ^ 0 then the change of variables of Equation (21.13) with C = P c changes the 
parametric curve to s 


X{t) = e Tlt 
Y(t ) = t Tr>t 


which leads to the implicit equation 


F(X, Y) = A r3 - y r ' = 0 


(21.14) 


with X , Y > 0 and t — -".W = 

~ ri r 3 * 


If 1^1 — 0 then P 0 is on one of the lines passing through P c in the direction of the 
eigenvectors of A, and either E x = 0 or E 2 = 0 (or Ei = E 2 = 0 if P 0 = P c ). The 

tangent curve will be a straight line. If E 2 = 0 then an implicit equation for the tangent 
curve is 


F(x, y) = (x - x c )e 2 i - (y - y c )e u = 0 

and in the case that E\ — 0 , we have 

F(x,y) = {x - x c )e 22 - (y - y c )e u - 0 

Note that these last two implicit equations are in terms of the original variables (a: u) See 
Figure 21.8. 
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Case 2. (0 = ri, r 2 ^ 0) 

If \E\ / 0, then the change of variables given by Equation (21.13) with C = Pq will 
yield 

X{t) = i 

y(t) = e tr2 - 1 
and 

F{X,Y) z=Y-e r2X + 1 = 0 (21.15) 


withX >0, Y > — 1, f = X. 

If \E\ — 0, then either 2? 2 = 0 and (in the original variables) we have the implicit 
equation 

( x — x 0 )e 2 i - (y — 2/o ) e 1 1 = 0 , 


or Ei = 0 and we have the equation 

(x - x 0 )e 2 2 - {y - 2/o)^i2 = 0 . 


See Figure 21.9. 
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Case 3. (rq = 0 = r 2 ) 

If | S'! 7^ 0, then Equation (21.13) with C = P 0 , can be used to obtain 

y - X 2 = 0. x > 0, Y > 0; t = X = s/Y. ( 21.1 

If |£*| = 0, then E-> = 0 and (in the original variables) we have the implicit equation 
(x - x 0 )e 2 i - (y - J/o)en = 0. 


Case 4. (n = r 2 ^ 0 ) 

If |£| ^ 0. *en the change of variables of Equation (21.13) leads to 

Xln(X) - rY = 0, X, Y > 0. t = - 

~ ‘ X 

If |£'| = 0, then the implicit equation (in the original coordinates) is 
(x - X c )(y 0 - y c ) - (y - y c )(x 0 - x c ) = 0. 


Case 5. (y q- A i, y — A i, A 0) 

The change of variable given by Equation (21.13) leads to 


y 

— = tan 
X 



/n(.Y 2 + F 2 ) 


(21.17) 


(21.18) 
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21.3 Incremental Methods 

The basic idea behind incremental methods consists of starting at the entry point P(0) = 
Pq and then successively computing P{U) 9 i — 1,2,3,... for increasing values of the 
parameter t and stopping when P(ti) leaves the triangle. These successive values that 
are computed along the way to finding the exit point can also be used as the basis for 
displaying the curve. Often it is efficient and convenient to use equally spaced values of 
the parameter; that is, £, = zAt, i = 1,2,3,.... At least this can be the default strategy 
and the value of At can be adjusted every so often, depending on accuracy estimates or the 
spacing requirements between P(t) and P(t + At). 

In this section, we cover two types of incremental methods. In Section 21.3.1 we dis- 
cuss briefly the application of conventional methods such as Euler s and Runge-Kutta to 
the special case at hand of linearly varying vector fields. In Section 2 1.3.2 we discuss some 
incremental methods for computing values on the curve which are based upon the exact so- 
lutions which were explicitly given in Section 21.2. Before we proceed, we wish to briefly 
discuss two topics that come up in the implementation of incremental methods. Namely, 
how to choose A t and how to test when the curve leaves the triangle. We take up the latter 
topic first. 

The functions fi(x, y), fj(x> y), and f k (x, y) of Equation (2 1.4) were carefully defined 
so that f a (x 9 y) > 0 for points (x, y) on the same side of as P a , a = ij , or k. This 
assumes that the points P„ Pj, and P k are listed in counterclockwise order. So by testing 
these functions, we can determine if the next point leaves the triangle and also which edge 
it leaves from and this can then be used to determine the data for the next triangle it enters. 
A particularly efficient way to compute all three values is based upon the identity 

/(x, y) = *(/(!, 0) - /(0,0)) + y(/(0, 1) - /(0,0)) + /(0,0) (21.19) 


where 

/ fi(x,y) \ 

/(*.y)=l fj( x <y) ■ 

\ fk{x,y) ) 

As we previously mentioned, it is sometimes useful to be able to approximately control 
the distance between P(t) and P{t+ At) by choosing an appropriate A t. Say, for example, 
it is desired that 

\\P(t + At)-P(t)\\K6. 


The mean value theorem implies 

P(t + At) - P{t) = AtP'(r), t < r < t + At. 


Equation (2 1 .6) then yields 

P(t + At) - P(t) = A t(A(P{r) + B), 


Now a local estimate as in 


S 

PP(<) + B)|| 


t < t < t + At. 


At = 
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can be used or an overall estimate as in 


maxmW, 11 ^ 11 , 11 ^ 11 } ' 

This last estimate is based upon the property that ||.A.P + B\\ is bounded by (max f 1 1 14 1 1 
1 1 Vj 1 1> 1 1 14 1|}) for P in the triangle of interest. 

21 .3.1 Conventional Methods Applied to Linearly Varying 
Vector Fields 

It is possible to treat Equation (2 1 . 1 ) as an ordinary differential equation and to use standard 
and conventional numerical methods to compute approximations to P(t). The most popular 
numerical methods in this context are incremental methods. The particular computational 
formulas for a sampling of these methods are: 

Euler’s: 


P(t + At) a P(t) + At ■ V(p{t)) 

Runge-Kutta 2 nd Order: 

P(t + At) a P(t) + i(Vi + V 2 ) 

Vi = At ■ V(P(t)) 

V 2 = At ■ V(P(t)+Vi ) 

P(t + At) = P(t) + 5(14 + At • V(P(t) + 14)) 

Runge-Kutta 4 £/1 Order: 

P(t + At) a P(t) + ±(14 + 2 Vj + 2 V 3 + V 4 ) 

Vi = At • V(p(t)) 

V2 = At • V(p(t) + ±V 0 
V 3 = At ■ V(p(t) + ±V' 2 ) 

V 4 = At • V(p{t) + ^V 3 ) 


Adaptive Runge-Kutta-Fehlberg: 


P{t + At) a P(t) + ^14 


+ 01^ 


Vi = At • V(p{t)) 


V 2 = At - V(p(t) + ± 14 ) 
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V 3 = At-V(p(t) + j i V 1 + ^V 2 ) 

V*=A t • V(p(t) + AfffVi - 

V5 = At • V(p(t) + - 8^ 2 + ¥§V 3 - f^V 4 ) 

V 6 = At- V(p(t) - AVi + 2 V 2 - ||gVs + if§| K, - &Vi) 

The above formulas can be simplified somewhat when the special property that V(P) = 
AP + B is utilized. The formulas in this case become: 


Euler’s: 

P(t + At) S (/ + A tA)P(t) + AtB 


Runge-Kutta 2 nd Order: 

P(t + At) = (/ + At A + i ^ 2 )P(t) + At(I+ ^ )B 


Runge-Kutta 4 t/l Order: 

P(t + At) = (E-= 0 /»(*) + At (£?=o P(*) 


Adaptive Runge-Kutta-Fehlberg: 


P(t + At) 


ao ( A t A ) 1 . 

l^iz rO i\ “ r 


+ A< (EL. 


(AM)' 

(i+i): 


(At-A ) 5 \ 

2080 y 


5 


Error Estimate = (At/l) s (^g - 7^) 


21 .3.2 Incremental Methods for the Exact Solution 


Here we are looking for a formula of the form 

P(t + At) = G{At)P(t) + H(At) 

which will allow computation of points at equal parameter values on the curve to be com- 
puted with only the application of an affine transformation. Repeated application of the 
basic differential equation which characterizes P(t ) leads to 


P'(t) = AP(t) + B 

P"(t) = AP'{t) = A 2 P{t) + AB 

P'"(t) = AP"{t) = A 3 P(t) + A 2 B 
which leads to the expansion 

P(t + At) = (I + AtA+ ( -^-+ -)P(t) 

+(At + (Af) 2 ^ + (Af) 3 ||- + - 
= E(At)P(t) + F(At)B 


)B 
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where E{At) = / + F(At)A. Both of these series for E(At) and F(At) have limits that 
can be approximately calculated in a variety of ways. One, of course, is to simply truncate 
the series. If we utilize the expressions for the solution given in Section 21.2.1, we arrive 
at the following particularly efficient methods which rely only upon the need to compute 
the functions e‘, sin(t), and cos(f). 

The matrix E referred to below is {E y , E 2 ) where these two vectors are defined in the 
different cases covered in Sections 2 1 .2. 1 and 2 1 .2.2. 


Case 1 . ( 0 ^ ri ^ r 2 ^ 0) 

If |JF| ^ 0 then 

p(t + At)-p c = E ( e ^ ri eA ° r , )fr-‘(P(o-p c) . 

If 1^1 = 0 ^en P 0 is on one of the lines passing through P c in the direction of the 
eigenvectors of A and either E x = 0 or E 2 = 0 (or E x = E? = 0 if P 0 = P c ). The tangent 
curve will be a straight line and so 

P(t+ Af) = P(t)±At(P e - P 0 ). 


Case 2. (0 = r it r 2 ^ 0) 


p(t + At) = [/ + { ±!2-Jl A] p {t) + At[I + 


r 2 


Atrl 


Case 3. (r, = 0 = r 2 ) 


P(t + At) = (I + AtA)P(t) + (Af + — )B. 


Case 4. {r x = r 2 ^ 0) 

If 1 2E7J ^ 0, then 

P(t + Af) -P c = Ee Atr ' ( ^ J ) E~HP(t) - Pc)- 

If (Fj = 0, then P(t) is a straight line and so 

P{t 4- At) = P(t) ± A t{P c - P 0 ). 

Case 5. (/j + A i, ji — A i, A ^ 0) 

Pit + At) - P, = Ee^ ( C0S ( AA< ) \ p) 

V 1 c { sin(XAt) cos(XAt) J b (P(t) ~ Pc)- 
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21 .4 Working Directly with Barycentric Coordinates 


In this section, we discuss the representation and computation of tangent curves using 
barycentric coordinates rather than conventional Cartesian coordinates. There are some 
potential advantages to the use of barycentric coordinates depending upon the particular 
application. One advantage to barycentric coordinates is that it is trivial to determine 
whether or not a point on a tangent curve is inside the triangular domain or not, as all 
three barycentric coordinates must be positive in this case. Another less obvious advan- 
tage is for applications where the triangular domain is not located in the normal xy-plane. 
We will show some examples of this later where the triangular domains consist of a trian- 
gulated approximation of the sphere. Here, even though the flow vectors are 3D vectors, 
they are in the plane of the domain and so all that we have developed in this chapter still 
applies. It would be possible to use a 3D affine transformation to move the problem to the 
plane, solve it, and transform back. But there is no need to do this transformation if we do 
all calculations in barycentric coordinates. But it should also be noted that most graphics 
routines require Cartesian coordinates and so the particular application environment has to 
be taken into consideration. 

For completeness, we give some background on barycentric coordinates. Given a point 



barycentric coordinates, b t , bj , and b k of this point relative to the triangle Tjjk with vertices 
P { , Pj , and Pk are defined by the relationships 


x 

y 


biPi + bjPj + bkPk 


( 21 . 20 ) 


1 = bj + bj + bk 


There are several alternative ways of defining or determining the barycentric coordinates 
of a point. For example, 



where A it A jt and A k represent the areas of the subtriangle shown in Figure 21.10 and A 
is the area of Tijk- Also, 


x — Xk 
y-yk 

Xj Xk 

Vj - Vk 

k . — _ 

X — Xi Xi — x k 

y - yi yi - yk 

Xj - Xk 

Xj - x k 

’> °J ~ 

Xj - X k Xi - x k 

Vi “ Vk 

yj - Vk 


1 

& 

1 

iS 


b k 


X — Xj Xj — Xj 

y - Vj Vi - yj 

X k -Xj Xj - Xj 

Vk - yj yi - yj 


( 21 . 21 ) 




Figure 21.10: Areas leading to barycentric coordinates. 


We can note that b u bj, and b k are simply linear functions of x and y. In fact, 6, is 
just the lmear function which is one at P { and zero at P 3 and P k . It will be convenient to 
explicitly denote this dependence so that we have 



and the relationship 



Also we will denote the column vector consisting of the three barycentric coordinates as 



It can be noted that 

6 ( 1 ) =X6 (o ) +y6 ( l l ) (21.22) 
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and 


b(aP + (1 - a)Q) = ab(P) + (1 - a)b{Q). (21.23) 

The 3x3 matrix whose columns consist of the barycentric coordinates of the three flow 
vectors at each of the vertices will be denoted by 


V = (b(Vi),b(Vjhf>(V k )). 


It should be noted that if 


/ ft 

Q = ft 
V 9fc 


has the property that <?,- +qj + q k = 1 (as will be the case when Q represents the barycentric 
coordinates of a point) then VQ will have this same property. It is also true that if g, + ft + 
q k - 0 then the entries of V Q will also sum to zero. 

It is easy to verify that 


A = {Pi,Pj,Pk)V[b 


-b 


(!)-‘(o)) (2L24) 


and 


A[ X y j + B = (Pi,Pj,P k )Vb[ X y 


(21.25) 


Also 

Vb {l)= ‘ |B > 

and if there exists a point P c such that 

AP C + B 


_ f 0 \ 
" V ° ) 


then 


6 ( 0 ) = vbm - 


We should point out that the 3 x 3 matrix V has the same eigenvalues as .4 and if AE = rE 
for a 2D vector E and scalar r, then 


V(b(E) - b 


= r(b(E)-b[ 


)• 


In addition V has the eigenvalue 1 . 
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The tangent curve P(t) which was previously given in Cartesian coordinates will be 
denoted in barycentric coordinates by 


/ MO \ 

p(0 = *(P(0)= WO 
\ MO J 

and the basic differential equation which characterizes P(t ) now takes the form 

P'(0 = Vp{t) ~ b ^ ® j (21.26) 


withp(O) = b{P 0 ). 

Using the same notation as in Equation (21.7) we can now represent the tangent curve 
which is the solution of Equation (2 1 .26) in barycentric coordinates by 


P(0 = <h(0 


b(Ei)-b 


0 

0 


+ $2(0 


b(E 2 ) - b 


0 

0 


+ b(C) (21.27) 


where, as before, each of the cases selected by the eigenvalues of ,4 (or V now) determine 
the bases functions 4>i(f) and <f- 2 (() and the vectors E u E 2 , and C. We now discuss 
incremental methods in terms of barycentric coordinates. In this context we are looking for 
a formula of the form 


p(t + At) = g(At)p(t) + h{At). (21.28) 

Similar to before, we can use repeated applications of Equation (21.26) to yield 

P(t + At) = (I + AtV+ ( -^pl + ..-)p(t) 

= + f{At)b ^ ° Q j . (21.29) 


Also, we can note that e(A t) = I -f- f(At)V. These series for e(A2) and f(At) can 
be approximated in a variety of ways. As we saw earlier, if we truncate the series, then 
we have the equivalent of an Euler or Runge-Kutta method applied to the special case of 
linearly varying vector fields. The same type of results holds in the case of barycentric 
coordinates. The results are: 
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Method 

e(At) 

/(At) 

Euler 

(I + AtV) 

—(At/) 

Runge-Kutta 
2 nd Order 

(TTAtv + fKFv 7 ) 

-(At/+ TAFV) 




Adaptive Runge- 
Kutta-Fehlberg 

(/ + AtV+ ±At'V* 

+ §At 3 l/ 3 + ±At 4 V 4 
+ihAt*V* + i&- 0 A t sV 6 ) 

-{Atl+^At-V 
+ |At 3 V r2 + ^At 4 V 3 
+ ^At*V 4 + ^At 6 V 5 ) 

Adaptive Runge- 
Kutta-Fehlberg 
error estimation 

(f^At^ + ^A t 6 V 6 ) 

-(^At 5 ^ + ^At 6 l/ 6 ) 


Similar to what was done in Section 21.3.2, we can develop exact representations uti- 
lizing transcendental functions. The details for the various cases follow. The matrix D 
referred to below is 

and D + = {D t D)~ 1 D t . The two vectors E\ and E 2 are defined in the different cases 
covered in Sections 2 1 .2. 1 and 2 1 .2.2. 

Case 1. (0 ^ n ^ r 2 ^ 0) 

Ifl^l ^ 0, then 

/ p Atr i rj \ 

p(t + At) - b(P e ) = D [ Q e Atr 3 J D + (p(t)-b(P c )). 

If \E\ = 0, then P 0 is on one of the lines passing through P c in the direction of the 
eigenvectors of A and either E\ = 0 or E 2 — 0 (or E\ = E 2 = 0 if Po = P c )- The tangent 
curve will be a straight line and so 

p(t + At) = p(t) ± At(b(P c ) - b(P 0 )). 


Case 2. (0 = ri,r 2 ^0) 



p(t + At) 
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CdS6 3. (Vi 0 — Vo) 

p(t + At ) = (f + AtV+(e* c -At- l)V 2 )p(t) 

-(Atl + — ^ - + (e At - -At- l)V 2 )b ^ jj j 

Case 4. (n = r 2 ^ 0) 

If \E\ ^ 0, then 

Pit + At) - HP.) = De ( ^ ( J ) C+WO - 

If 1^1 = 0, then the tangent curve is a straight line and so 

p(t + A*) = p(t) ± A t(b(P c ) - b(P Q )). 


CdS6 5, (p + A i, p — A i, A ^ 0) 
p(t + AO - b(P c ) 

_ / cos(AA^) -sm(AA^) \ _ , . , % 

\ sin(XAt) cos(XAt) 1 ) D+ (p(t)-b(P c )). 

21.5 Topological Methods 

Topological methods provide a clear and uncluttered means of visualizing a two- 
dimensional vector field. They give a good overview of the flow with relatively little in- 
formation, but they can require some effort to compute. They were first introduced into 
the visualization literature by Helman and Hesselink [17]. They consist of a collection of 
tangent curves which separate the flow into regions. The tangent curve boundaries of these 
regions connect together certain critical points. Critical points are locations where the flow 
is zero. 

In a nutshell, there are two main steps to computing a topological graph. First, all 
critical points are computed and classified on the basis of the nature of the flow near the 
critical point. In the second step, certain tangent curves are started at critical points and 
traversed out until they either link up with other critical points or leave the domain. In this 
section, we will discuss what is necessary to compute a topological graph within the special 
context of this chapter which is piecewise linear vector fields over triangulated domains. 
An example of a topological graph is shown in Figure 21.11. 

We now discuss some details of the critical point computation and classification step of 
the algorithm. In a general context with the vector field 


V{x,y) 


_ ( u(x,y) 

v(x,y) 


a critical point, P c , is simply a position where 

u (x c , Vc) 

v(x cVc) 


V(Pc) = 


0 

0 



21 .5 Topological Methods 


541 



Figure 21.11: A topological graph of a two dimensional vector field. 


The local behavior close to a critical point is determined by the Jacobian 

U \ _ ( u A x c,Vc) «y {Xc,y C ) \ 

(£e,yd ^ v x (x c ,y c ) v y (x c ,y c ) )' 

This is due to the fact that the flow field is approximated by the linear terms of its expansion 
near the critical point. That is, 

V(P ) 2 V(P C ) + {P- Pc)J{Pc) = {P- Pc)J(Pc). 

Using techniques similar to those used in Section 2 1 .2. 1 , it is determined that there are six 

different types of critical points (see Figure 21.12): .... 

i) Saddle Point J {Pc) has two real, nonzero eigenvalues which differ in sign 

ii) Attracting Node J{P C ) has two negative eigenvalues 

iii) Repelling Node J(P e ) has two positive eigenvalues 

iv) Attracting Focus J(P C ) has complex eigenvalues /i + Ai, /r - Ai, A ^ 0, /r < 0 

v) Repelling Focus J{P e ) has complex eigenvalues y + A i, y. - \i, A / 0, y > 0 

vi) Center J(P C ) has purely imaginary eigenvalues +Ai and -A i, A ^ 0 

In general, the computation of critical points can be a rather formidable problem. For 
example, in the case of curvilinear grids, cells must be searched for possible candidate 
cells and then a nonlinear system (two equations in two unknowns) of equations must 
be solved. This normally requires some iterative method which could possibly fail to 
converge or converge to a point which is actually not a critical point. In the case of a 
linearly varying vector field, the situation is much simpler. Only the linear system 



Xc 

Vc 


+ 5 = 0 
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Figure 21.12: Characterization of critical values. 


must be solved. Note that .4 and B are defined the same as earlier in Equation (21.2). We 
next determine if this is a “real” critical point by deciding whether or not P c is actually 
in the triangle. This determination can be made by computing the barycentric coordinates 
of Pc (or possibly a scaling of them) and checking to see if all three of them are positive. 
More discussion on barycentric coordinates can be found in the previous Section 21.4. 

Often, a flow field will have an interior boundary which represents an object about 
which the flow is being analyzed. All the points on this interior boundary are critical points 
as the flow is forced to be zero here. Some of these points are of special interest. These are 
the points of attachment (at) and detachment (de) which belong to tangent curves which 
separate the flow along and near the inner boundary. In some previous discussion in the 
literature, the characterization of these points has not been so rigorous, but here in the 
context of piecewise linear flow fields, we can be very precise. Consider a triangle with 
only one vertex on the inner boundary. This vertex will necessarily be a critical point. If the 
eigenvalues of the Jacobian for this triangle indicate that this critical point is classified as a 
saddle point, then it is a candidate for a point of attachment or detachment. It will be a point 
of attachment if the eigenvector (or its negative) associated with the negative eigenvalue 
lies between the two edges of the triangle sharing this critical point. If an eigenvector 
associated with the positive root lies in the triangle, then this critical point will be a point 
of detachment. See Figure 21.14. 

It should be noted that it is possible for a single point to be both a point of attachment 
and a point of detachment. This is illustrated in Figure 21.15. We should also point out that 
it is possible for a single point on the inner boundary to belong to two different triangles 
which have only this point on the inner boundary and because of this, this single point 
could be classified as a certain type of critical point for one triangle and another (or the 
same) type of critical point as a vertex for a different triangle. An example of this is shown 
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in Figure 21.16 where the critical point on the upper-left portion of the inner boundary is a 
point of attachment for one triangle and a point of detachment for an adjacent triangle. 

The second major step in computing a topological graph consists of the linking al- 
gorithm. From each saddle point which is interior to a triangle* four tangent curves will 
emanate — two associated with each ei genvector of the Jacobian. One curve emanates in the 
direction of the eigenvector and one in the negative direction of the eigenvector. The curves 
emanating in the directions of the eigenvector associated with the positive eigenvalue will 
be traversed in positive parameter direction and these curves will move along with the flow. 
They have the chance to link up with attracting nodes or foci or possibly other saddle points 
along the eigenvectors of inward flow to the saddle point. The two curves emanating in the 
direction of the eigenvectors associated with the negative eigenvalue will be traversed in 
a negative parameter direction and will move along in a direction opposite (or negative) 
to the direction of the flow field. They have the chance to link up with repelling nodes or 
foci or other saddle points. From each point of detachment, one curve will emanate and be 
traversed in positive parameter direction. From each point of attachment, one curve will 
emanate and be traversed in a negative or opposite direction to the flow. The algorithm is 
complete when each of these curves has linked to another critical point or leaves the domain 
either by encountering the outer boundary or the inner boundary. An example of the results 
of this algorithm are shown in Figure 2 1 . 1 6. We should point out that this linking algorithm 
does not produce all separating tangent curves for it is possible that in addition to the center 
points, some repelling or attracting nodes or foci could be left unconnected to any tangent 
curve. However, this does not occur when the domain is a triangulated approximation to 

the sphere. 
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21.6 Examples 

The first example that we include uses a simple quadratic equation to specify the vector 
field. The equations are given in Equation (21.30). One reason for including this type of 
example is to allow other implementors to easily verify and compare their software results. 


U{x, y) 


V{x,y) 


-0.103209 + 0.05151 lx - 0.302699;/ 

+0.037546x2/ — 0.232875x 2 + 0.61 1528r/ 2 

(21.30) 

0.143656 + 0.687847x - 0.144779;/ 

-0.213010x2/ - 1 029676x 2 + 0.246278;/ 2 


As we have mentioned earlier, the methods described here can be applied to any tri- 
angulated domain. In Figure 21.18 we show the topological graph and some additional 
tangent curves (in cyan) for a vector field defined over a triangulation of a spherical do- 
main. Similar to Figure 21.17, two different resolutions of the triangulation are shown. In 
the right column, each triangle of the left column has been replaced by four subtriangles. 
We have intentionally used flat shading rather than Gourard shading for the rendering of 
the sphere so that the triangulation is apparent. 

In the next example, we illustrate the use of the methods developed here to visualize 
a multiresolution model of a vector field over a curvilinear grid. The topological graph is 
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Figure 21.15: A point on the inner boundary which is both a point of attachment and a 
point of detachment. 


shown at four different levels of approximation in Figure 21.19. One unique feature allowed 
by these methods is that no matter how coarse the resolution becomes, the boundaries 
have not changed from the original data. This is a particularly important aspect for the 
inner boundary which is often the focus of attention for a flow analysis. In Figure 2 1 .20, 
we show the topological graph for some partial reconstructions of the flow field. These 
methods allow the user to zoom in and out and only reconstruct the portion of interest. 

In Figure 21.21 we show results similar to those of Figures 21.19 and 21.20 except that 
now the domain is a triangulated sphere and the reconstruction is not done on the basis 
of regions, but on the basis of the magnitude of the coefficients of the wavelet basis func- 
tions. This data represents “real” data provided to us by Roger Crawfis and Nelson Max 
of Lawrence Livermore National Laboratory. It is one time step and one tier of simulated 

wind velocity data. 

The data used for the examples shown in the images of Figure 21.22 and Figure 21.23 
was provided to us by Yasuo Nakajima, Nissan Motor Company, Japan. Actually, the 
domain of this data is three-dimensional and not two-dimensional as required for the algo- 
rithms covered here. Figure 21.23 shows one slice through this 3D data. The 3D velocity 
vectors were projected into the plane of the slice. Many of the results of this chapter have 
been extended to 3D domains where the vector field is assumed to vary linearly over a tetra- 
hedrization of the domain. See [33] for more discussion on tetrahedrizing a 3D curvilinear 
grid. The tangent curves and critical points shown in Figure 21.22 were computed using 
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Figure 21.16: Topological graph for a curvilinear grid. A point on the inner boundary is 
both a point of attachment and a point of detachment. 
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TRIANGLE METHOD 5x5 TRIANGLE METHOD 15x15 



Figure 21.17: Explicit method used to compute topological graph for vector field given by 
Equation (2 1 .30). Two different resolutions for the triangulation are shown. 

extensions of the algorithms discussed here. See [22] for more details. 


21.7 Conclusions and Remarks 

1 . One of the advantages of the explicit methods we have developed here is that accu- 
racy can be monitored and controlled. Runge-Kutta and other incremental methods 
for solving ODEs are notorious for “wandering off” the true solution and once an 
error is made, there is no way to recover because from the erroneous point forward, 
the method is attempting to solve a different (wrong) problem than the one it set out 
to solve in the beginning. Variable step size methods which estimate the error and 
attempt to control it by adaptively changing the step size are helpful in this regard, 
but the problem is that the error is only estimated and only guesses about the prox- 
imity of the computed solution to the desired solution can be made. This is the case 
for the R-F-K method covered in Sections 21.3.1 and 21.4. In the proposed method, 
the accumulated error is a result of how accurately the intersection of a cell boundary 
and a particular explicitly defined tangent curve are computed. This is formulated as 
a root finding problem and so this computation can be done as accurately as deemed 
necessary. 

2. Another advantage of the present method is speed and efficiency. Since the tangent 
~ curve is known explicitly for each triangle, parameters pertaining to this definition 

can be precomputed and stored and then used when a particular tangent curve pene- 
trates this triangle. The global shape of the tangent curve is known by its sequence 
of entry and exit points for each triangle it intersects and so the overall appearance 
of the curve is not seriously degraded if a local approximation for each cell is used. 
We have used a parametric cubic Hermite curve on each triangle with good success. 
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Figure 21.18: Explicit method used to compute critical values and tangent curves for a 
vector field defined over a spherical domain. 
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Figure 21.21: Topological graphs for wind data over the earth are computed using the 
explicit method. The right column is wavelet reconstruction based on the largest 3% of the 
wavelet coefficients. 
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Figure 21.22: The explicit method is used to compute tangent curves linking critical values 
for a 3D vector field. 
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3. Because of the general nature of the approach of the methods covered here, they can 
be extended to any domain which consists of a collection of triangles. This includes 
manifolds of arbitrary topological type if a triangulated approximation of the domain 
is acceptable. 


4. The computation of critical points is greatly simplified. Critical points are crucial 
to topological methods. If the vector field which has been defined by a particular 
cell and extended to the entire domain of E 2 has a critical value, then it can be 
computed by solving a linear system of equations. A solution to the linear system 
of equations is a “rear critical value for the piecewise defined vector field if it lies 
within the cell, otherwise it is not. Normally, the computation of critical values 
requires a rather tedious computation involving testing whether or not a cell might 
have a critical value, followed by the solution to a nonlinear system of equations by 
Newton’s methods or some other iterative scheme. See [17], for example. 


5. The last advantage we mentioned is that the present method does all of its compu- 
tations in physical coordinates as opposed to computational coordinates. To some 
this might initially appear to be a disadvantage since computational coordinates are 
introduced for the specific purpose of their namesake. But if the domain is triangu- 
lated and linear variation is assumed over each triangle, we not only gain a method 
that is affine invariant, but there is no need to map the data to computational coordi- 
nates, solve the problem and map the solution back; we simply compute the solution 
directly in physical space. 


6. We first mentioned the basic ideas of these methods and reported on some prelimi- 
nary results at a presentation given in 199 1 at the first Dagstuhl Seminar on Scientific 
Visualization. At this time, Nelson Max pointed out that he had mentioned the idea 
in an earlier paper [29]. Sawada [42] has also mentioned the idea of using explicit 
representations. We presented some preliminary results on the 3D extension of the 
methods covered here at a technology assessment workshop held the Summer of 
1993 in Darmstadt, Germany. See [40] for the proceedings of this workshop. Fig- 
ure 21.22 has previously appeared in [34]. 


Acknowledgments 

This research was supported by the National Aeronautical and Space Administration under 
NASA-Ames Grant, NAG 2-990. Partial support was also provided by the North Atlantic 
Treaty Organization under grant RG0097/88. The data for the example of Figures 21.22 
and 21.23 was provided by Yasuo Nakajima of Nissan, Japan. The data for the example 
of Figure 21.21 was provided by Roger Crawtis and Nelson Max of Lawrence Livermore 
National Laboratory. 



BIBLIOGRAPHY 


555 


Bibliography 

[1] F.H. Bertrand and P.A. Tanguy, “Graphical Representation of Two-Dimensional 
Fluid Flow by Stream Vectors,” Communications in Applied Numerical Methods, 
Vol.4, 1988, pp. 213-217. 

[2] J. Buckmaster, “Perturbation Technique for the Study of Three-Dimensional Sepa- 
ration,” Phys. Fluids, Vol. 15, 1972, pp. 2106—21 13. 

[3] P.G. Buning, “Numerical Algorithms in CFD Post-Processing,” von Karman Insti- 
tute for Fluid Dynamics, Lecture Series 1989-07. 

[4] B.J. Cantwell, “Similarity Transformations of the Two-Dimensional Unsteady 
Stream Function Equations,” / Fluid Mech., Vol. 85, 1978, pp. 257-271. 

[5] M.S. Chong, A.E. Perry and B.J. Cantwell, “A General Classification of Three- 
Dimensional Flow Fields,” Physics of Fluids A, Vol. 2, No. 5, 1990, pp. 765-777. 

[6] D. Darmofal and R. Haimes, “An Analysis of 3D Particle Path Integration Algo- 
rithms,” Proc. 1995 AIAA CFD Meeting, San Diego, Calif., 1995. 

[7] A. Davey, “Boundary Layer Flow at a Point of Attachment,” / Fluid Mech., Vol. 10, 
1961, pp. 593-610. 

[8] R. Denzer, “Application of Visualization in Environmental Protection,” Focus on 
Scientific Visualization, H. Hagen, M. Muller and G.M. Nielson, editors, Springer, 
1993, pp. 73-82. 

[9] R.R. Dickinsion, “Interactive Analysis of the Topology of 4D Vector Fields,” IBM 
Journal of Research and Development, Vol. 35, No. 1/2, 1991, pp. 59—66. 

[10] D.S. Ebert and R.E. Parent, “Rendering and Animation of Gaseous Phenomena by 
Combining Fast Volume and Scanline A-Buffer Techniques, Computer Graphics, 
Vol. 24, No. 4, 1990, pp. 357-366. 

[11] M. Fruehauf, “Combining Volume Rendering with Line and Surface Rendering, 
Eurographics ’91, F.H. Post and W. Barth, editors. North Holland, 1991, pp. 21-32. 

[ 1 2] R.S. Gallagher, “Span Filtering: An Optimization Scheme for Volume Visualization 
of Large Finite Element Models,” Proceedings of Visualization ’91, G.M. Nielson 
and L. Rosenblum, editors, IEEE Computer Society Press, 1991, pp. 68-75. 

[13] R.S. Gallagher and J.C. Nagtegaal, “An Efficient 3-D Visualization Technique for 
Finite Element Models and Other Coarse Volumes,” Computer Graphics, Vol. 23, 
No. 3, 1989, pp. 185-194. 

[14] M. Geiben and M. Rumpf, “Visualization of Finite Elements and Tools for Numeri- 
cal Analysis,” Advances in Scientific Visualization, F.H. Post and A.J.S. Hin, editors, 
Springer, 1992. 


556 


Tools for Computing Tangent Curves and Topological Graphs 


[15] R.B. Haber and D.A. McNabb, “Visualization Idioms: A Conceptual Model for Sci- 
entific Visualization Systems,” Visualization in Scientific Computing, G.M. Nielson, 

L. Rosenblum and B. Shriver, editors, IEEE Computer Society Press 1990 pp 74— 
92. 

[16] B. Hamann, D. Wu, and R. Moorhead, “Flow Visualization with Surface Particles,” 
IEEE Transactions on Visualization and Computer Graphics, Vol 1 No 3 Sep 

1 995, pp. 210-217. ' 

[17] J Helman, and L. Hesselink, “Representation and Display of Vector Field Topology 
in Fluid Flow Data Sets,” IEEE Computer, Vol. 22, No. 8, 1989, pp. 27-36. 

[18] J. Helman and L. Hesselink, “Surface Representation of Two- and Three- 
Dimensional Fluid Flow Topology,” Proceedings Visualization '90, A. Kaufman, 
editor, IEEE Computer Society Press, 1990, pp. 6-13. 

[19] J. Helman and L. Hesselink, “Visualizing Vector Field Topology in Fluid Flows,” 
IEEE Computer Graphics and Applications, Vol. 11, No. 3, pp. 36-46. 

[20] W. Hibbard and D. Santek, “Visualizing Large Data Sets in the Earth Sciences,” 
IEEE Computer, Vol. 22, No. 8, pp. 53—57. 

[21] S. Hultquist, “Interactive Numeric Flow Visualization Using Stream Surfaces,” 

Computer Systems in Engineering, Vol. 1, Nos. 2-4, 1990, pp. 349-353. 

[22] I.-H. Jung, Topological Visualizing Method of Vector Fields in Fluid Data Sets, M.S. 
thesis, Arizona State University, Department of Computer Science and Engineering 
June 1993. 

[23] J.C.R. Hunt, C.J. Abell, J.A. Peterka and H. Woo, “Kinematical Studies of the Flows 
Around Free or Surface-Mounted Obstacles; Applying Topology to Flow Visualiza- 
tion,”^ Fluid Mech., Vol. 86, 1978, pp. 179-200. 

[24] D.N. Kenwright, Dual Stream Function Methods for Generating Three-Dimensional 
Stream Lines, Ph.D. thesis, University of Auckland, Department of Mechanical En- 
gineering, Aug. 1993. 

[25] D.N. Kenwright and D.A. Lane, “Interactive Time-Dependent Particle Tracing Us- 
ing Tetrahedral Decompostion,” IEEE Transaction on Visualization and Computer 
Graphics, Vol. 2, No. 2, June 1996, pp. 120-129. 

[26] M.J. Lighthill, “Attachment and Separation in Three-Dimensional Flow,” Laminar 
Boundary Layers, L. Rosenhead, editor, Oxford University Press, 1963, pp. 72-82. 

[27] R. Lohner and J. Ambrosiano, “A Vectorized Particle Tracer for Unstructured Grids,” 

J. Computational Physics, Vol. 91, 1990, pp. 22—31. 

[28] K.-L. Ma and P.J. Smith, “Cloud Tracing in Convection-Diffusion Systems,” Pro- 
ceedings of Visualization '93, G. Nielson and L. Rosenblum, editors, IEEE Com- 
puter Society Press, 1993, pp. 253-260. 



BIBLIOGRAPHY 


557 


[29] N. Max, Private Communication , 1991. 

[30] W. Merzkirch, Flow Visualization , Academic Press, 1987. 

[31] C. Moler and C. Van Loan, “Nineteen Dubious Ways to Compute the Exponential 
of a Matrix,” SIAM Review, Vol 20, No. 4, Oct. 1978, pp. 801-836. 

[32] G.M. Nielson, “Modeling and Visualizing Volumetric and Surface- on-Surface Data ” 
Focus on Scientific Visualization, H. Hagen, M. Muller and G.M. Nielson, editors, 
Springer, 1993, pp. 191-240. 

[33] G.M. Nielson, “Tools for Triangulation and Tetrahedrization,” Chapter 20 in this 
volume. 

[34] G.M. Nielson and A. Kaufman, “Visualization Graduates,” Computer Graphics and 
Applications, Vol. 14, No. 5, Sep. 1994, pp. 17—18. 

[35] T.V. Pathomas, J.A. Schiavone and B. Julesz, “Applications of Computer Graphics 
to The Visualization of Meteorological Data,” Computer Graphics, Vol. 22, No. 4, 
pp. 327-335. 

[36] A.E. Perry and B.D. Fairlie, “Critical Points in Flow Patterns,” Adv. Geophys., Vol. 
18B, 1974, pp. 299-315. 

[37] A.E. Perry and D.K.M. Tan, “Simple Three-Dimensional Motions in Coflowing Jets 
and Wakes,” J. Fluid Mech., Vol. 141, 1984, pp. 197-231. 

[38] A.E. Perry and M.S. Chang, “A Description of Eddying Motions and Flow Patterns 
Using Critical-Point Concepts,” Ann. Rev. Fluid Mech., Vol. 19, 1987, pp. 125—155. 

[39] F.H. Post and T. van Walsum, “Fluid Flow Visualization,” Focus on Scientific Vi- 
sualization, H. Hagen, M. Muller and G.M. Nielson, editors. Springer, 1993, pp. 
1-40. 

[40] L. Rosenblum et al., editors. Scientific Visualization: Advances and Challenges, 
Academic Press and IEEE Computer Society Press, 1994. 

[41] A. Sadarjoen, T. van Walsum, A.J.S. Hin, and F.H. Post, “Particle Tracing Algo- 
rithms for Curvilinear Grids,” Proceedings Fifth Eurographics Workshop on Visual- 
ization in Scientific Computing, Rostock, Germany, May 1994. 

[42] K. Sawada, “Visualization of Unsteady Vortex Motion in the Flow over a Delta 
Wing,” Proc. of the 5th Int. Symp. on Computational Fluid Dynamics, Sendai, Vol. 
Ill, 1993, ISCFD. 

[43] W.J. Schroeder, C.R. Volpe and W.E. Lorensen, “The Stream Polygon: A Technique 
for 3D Vector Field Visualization,” Proceedings Visualization '91, L. Rosenblum 
and G.M. Nielson, editors, IEEE Computer Society Press, 1991, pp. 126-132. 

[44] N. Srinivasan, An Efficient Method for the Computation of Tangent Curves, M.S 
thesis, Arizona State University, Department of Computer Science and Engineenng, 
May 1995. 



558 


Tools for Computing Tangent Curves and Topological Graphs 


[45] J. Stolk and J.J. van Wijk, “Surface-Particles for 3D Flow Visualization,” Advances 
in Scientific Visualization, Springer, 1992. 

[46] T. Strid, A. Rizzi and J. Oppelstrup, “Development and Use of some Flow Visualiza- 
tion Algorithms, von Karman Institute for Fluid Dynamics, Lecture Series 1 989-07. 

[47] M. Tobak and D.J. Peake, “Topology of Two-Dimensional and Three-Dimensional 
Separated Flows,” AIAA 12th Fluid and Plasma Dynamics Conference Williams- 
burg, Va„ July 23-25, 1979. 

[48] M. Tobak and D.J. Peake, “Topology of Three-Dimensional Separated Flows” Ann. 
Rev. Fluid Mech., Vol. 14, 1982, pp. 61-85. 

[49] S.K. Ueng, K. Sikorski, and K.-L. Ma, “Efficient Streamline, Streamribbon, and 
Streamtube Constructions on Unstructured Grids,” IEEE Transactions on Visualiza- 
tion and Computer Graphics, Vol. 2, No. 2, 1996, pp. 100-1 10. 

[50] J.J. van Wijk, “A Raster Graphics Approach to Flow Visualization,” Eurographics 
'90, D.A. Duce and C.E. Vandoni, editors, North Holland, pp. 251-259. 

[51] J.J. van Wijk. “Spot Noise — Texture Synthesis for Data Visualization,” Computer 
Graphics , Vol. 25, No. 4, 1991, pp. 300-318. 

[52] W.J. Yang, editor. Handbook of Flow Visualizaton, Hemisphere Publishing, 1989. 



Chapter 20 

Tools for Triangulations and 
Tetrahedrizations and 
Constructing Functions 
Defined over Them 


Gregory M. Nielson 


20.1 Introduction 

This chapter is about triangulations and tetrahedrizations and functions defined over them. 
The original and main goal was to provide some information about tetrahedra and tetra- 
hedrizations and functions defined over them, but it was quickly realized that many of 
these topics are easier to describe and understand with some background on their two- 
dimensional analogs. Therefore, it was decided to also include material on triangulations. 
While some of the material exists elsewhere in the literature, much is new and appears here 
for the first time. The intended purpose for this chapter is to serve as a survey/tutorial in 
the area of data modeling and visualization. As data modeling and visualization becomes 
more sophisticated in its application domains and begins to deal with data sets that are more 
complex than Cartesian grids, there will be the need for tools to deal with these data sets. 
We feel that the tools and techniques covered here are very basic and will prove to be useful 
in a variety of contexts in data visualization. 

Now we have some comments about the organization of this chapter. While tetra- 
hedrizations are the goal, researchers have dealt with triangulations for a much longer pe- 
riod of time than tetrahedrizations and so triangulations and related matters are much better 
understood. The next section is a survey of triangulations and related matters of interest in 
modeling and visualization. The following section is on tetrahedrizations and we attempt 
to follow the same flow of information as in the section on triangulations as well as pos- 
sible. We use the phrase “as well as possible” because some aspects of triangulations do 
not generalize to tetrahedrization and some facts known about triangulations and triangu- 
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lar domains are yet to be known about tetrahedrizations and tetrahedral domains. On the 
other hand, there are topics of interest to tetrahedrization which have no 2D counterpart of 
interest — for example, visibility sorting for tetrahedrizations. The outline of this chapter is 
very simple. In Section 20.2 we go through a list of topics on triangulations and triangular 
domains and then in Section 20.3 we repeat these topics with reference to tetrahedrizations 
and tetrahedral domains. 


20.2 Triangulations 

20.2.1 Basics 

Definitions, Data Structures, and Formulas for Triangulations 

In order to avoid any possible confusion and problems later, it is usually best to be a little 
precise and formal about the definition of a triangulation. We start with a collection of 
points in the plane, P = (p, — ( Xi , t/,), i — 1, . . . , A r } and a domain of interest, D, which 
contains all of the points of P. We assume that the boundary of D is a simple (does not 
intersect itself), closed polygon. Often D is the convex hull of P, but in general, it need 
not be convex. In fact the boundary does not have to be a single polygon so that the 
domain is not even simply connected. ( Connected means that there is a path joining any two 
points and simply connected means that the complement is connected.) Roughly speaking, 
a triangulation is a decomposition of D into a collection of triangles which are formed from 
vertices of P. Since we are eventually interested in defining functions over D in a piecewise 
manner over each triangle, we must require that the triangles do not overlap so as not to 
have any ambiguities. Thus we require the collection of triangles of the triangulation to 
be mutually exclusive and collectively exhaustive. In order to continue this formalism to 
a precise definition, we need some additional notation. A single triangle with vertices p,-, 
Pj, and pk is denoted by Tijk and the list of triples which represents the triangulation is 
denoted by I t . A triangle T t j k is a closed 2D point set that includes its three edges which 
comprise its boundary. The interior of denoted by In^T^) is open and does not 
include the boundary. The edge joining p, and pj is denoted by and N e = {ij : ijk in 
I t for some k} is used to refer to the collection of all edges. Formally, the definition of a 
triangulation requires the following: 

i) No triangle T^k, ijk € h is degenerate. That is, if ijk £ I t then p it pj and p k are 
not collinear. 

ii) The interior of any two triangles do not intersect. That is, if ijk £ I t and aj3 7 £ I t 
then Int {T ijk )n Int(T ft/?7 ) = 0. 

iii) The boundary of two triangles can only intersect at a common edge. 

iv) The union of all the triangles is the domain D = Uijkeit^ijk- 

Examples of valid triangulations are shown in Figure 20.1 and Figure 20.2. Note that 
the example of Figure 20.1 is not convex and that of 20.2 is not simply connected. Even 
though the diagrams of Figure 20.3 and Figure 20.4 look all right, the actual triangulations 
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Figure 20.1: A triangulation of a nonconvex domain. 


given by the corresponding It ’s do not represent valid triangulations. In the case of Fig- 
ure 20.3 the triangle T 46 5 is degenerate. Even if this triangle is eliminated, what remains 
is not a valid triangulation because condition iii) would then be violated since edge e 4 6 
contains p*>. This example would become a valid triangulation if the point p$ were to be 
moved slightly to the right so as not to be on the edge e 4 6 - The information of Figure 20.4 
is not a valid triangulation because condition ii) is violated. 

We now want to make some assertions about the possibility of triangulating a domain 
containing a collection of data points that is bounded by a simple, closed polygon. First we 
note that in the case that the domain contains no interior data points, it is always possible 
to form a triangulation. Just for the sake of interest, we mention two ways that this can be 
accomplished. The first way is based upon the fact that every simple closed polygon with 
more than three vertices can be split into two polygons. This leads to an algorithm that 
recursively splits each subpolygon until only triangles are left. The following argument 
which guarantees that each simple closed polygon has a diagonal has been discussed in 
[16]. A diagonal is an edge between two vertices that lies inside the polygon and does not 
intersect the polygon except at the endpoints. 


Splitting a polygon: Let b be the vertex with minimum x-coordinate and ab and be be 

its two incident edges as is shown in Figure 20.5. If ac is not cut by the polygon, then uc is 
a diagonal. Otherwise there must be at least one polygon vertex inside T a bc- Let d be the 
vertex inside abc furthest from the line through a and c. Now edge bd cannot be cut by the 
polygon, since any edge intersecting bd must have one endpoint further from line ac.'The 
second approach leads to an iterative algorithm. We first give a definition. A vertex, p it of 
a simple, closed polygon is called protruding , provided the following conditions hold; 
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Figure 20.2: A triangulation of a domain which is not simply connected. 



1 4 2 

2 4 6 

2 6 3 

4 6 5 

4 7 5 

5 7 8 

5 8 6 


Figure 20.3: Not a valid triangulation. 
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Figure 20.6: Insertion of an interior point. 


i) The interior angle 0,-, between the edges, Pi-iPi and pipi+i is less than or equal to 
7r. (Cyclic notation is used here so thatp^v+i = p l .) 

ii) The triangle contains no other vertices of the polygon than pt-upi or 

P*+ i- 

iii) The interior of , l + l is contained in the interior of D. 

It is an easy matter to prove that every simple, closed polygon has at least one protruding 
vertex. (The proof is left to the reader. Some people call them ears and so there must be 
two of them!) We can triangulate the polygon-bounded domain by successively removing 
protruding vertices. This approach to triangulating the region bounded by a simple closed 
polygon is called the “boundary stripping algorithm.” It is easy to implement, but in a 
theoretical sense, it is not competitive with other algorithms (see, for example, the papers 
of Narkhede and Manocha [175] and Fournier and Montuno [94], among others). 

Once the boundary of D has been triangulated, it is a relatively simple matter to build 
a triangulation including the interior points. This can be done by simply inserting them 
sequentially in a manner which we now describe. 

Insertion of an interior point: If the point to be inserted, p, lies in the interior of the 
triangle T abc , we replace T abc with the three triangles: T abp , T bcp , T cap . Ifp lies on an edge 
shared by T abc and T bad , then we replace the two triangles T abc and T bad with the four 
triangles T bcp , T dbp ,T pca) T pad . These two cases are illustrated in Figure 20.6. 

It is also possible to generalize the insertion idea to include an edge. Once we are armed 
with this capability, we know that we can triangulate any polygon-bounded domain: simply 
connected or multiply connected (that is, with holes). 

Insertion of an interior edge: Assume that the one endpoint, p y lies in the triangle 
Tabc and that the other endpoint, q , lies in the triangle T xyz . See Figure 20.7. Collect all 
of the triangles from T abc to T xyz which are intersected by edge pq and form a region R 
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Figure 20.7: Insertion of an interior edge. 


with polygon boundary D. We can split D with polygon dpqw, where d is the vertex of T abc 
not on the edge common with the other triangles whose union is R, and w is the analogous 
vertex of T xyz . Now we know that each of these two domains can be triangulated. The 
union of these two triangulations, which each contain the edge pq, can replace the previous 
triangulation of D. 

In addition to I t , which represents the triangulation, it is often worthwhile to generate 
and maintain some auxiliary information about the neighbors of each triangle. This infor- 
mation is useful for traversal algorithms and evaluation algorithms which have a search- 
ing component that determines the particular triangle containing a point where a function 
defined piecewise over the triangulation is to be evaluated. One very common and par- 
ticularly useful data structure is that which is illustrated in Figure 20.8. The first three 
columns contain the data of /(, with the additional constraint that in reading from left to 
right (cyclically), the vertices of each triangle are traversed in a clockwise order. The next 
three columns contain the indices of the triangles which are neighbors to this triangle. The 
character <j> indicates that the triangle has an edge that is part of the boundary of D. The en- 
tries of these three columns are also in a special order. The fourth column contains the index 
of the triangle which shares the common edge with vertex indices specified in the second 
and third columns. Similar relationships hold for the fifth and sixth columns. The infor- 
mation represented by this data structure is called a “triangular grid.” The neighborhood 
information contained in the last three columns does not contain any “new” information 
over that of I t , but it is often the case (and this depends, of course, on the application) that 
it is useful data which is worth generating a priori. 

Another data structure for representing a triangulation which is useful for some ap- 
plications is illustrated by the example shown in Figure 20.9, which represents the same 
triangulation as that of Figure 20.8. Here, for each vertex, a list of all vertices which are 
joined by an edge of the triangulation is given. This list is given in counterclockwise order 
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Figure 20.8: An example that defines a triangular grid structure . 

Vertex Joining Vertices 

1 2, 3, 4, 5 

2 8, 7, 3, 1 

3 1,2, 7, 6, 4 

4 3, 6, 5, 1 

5 1,4, 6,7 

6 3,7, 5, 4 

7 6,3, 2, 8, 5 

8 2, 7 

Figure 20.9: The data that defines the data point contiguity list. 


around each vertex. This is called the data point contiguity list. We mention this particular 
data structure because of its convenience for dealing with the optimal Delaunay triangula- 
tion discussed in the next section. Also, it is very useful for computing the parameters of 
the Minimum Norm Network method [179], which is one of the most effective C 1 interpo- 
lation methods for scattered data. 

Even though there are a number of possible triangulations for any given domain D, the 
number of triangles is fixed once the boundary has been specified. More precisely, if N \ 
represents the number of vertices on the boundary and Ni the number of interior vertices 
so that N = Nt, + Ni, then the following formulas hold: 

N t = 2 Ni + N b -2 

and 


A r e = 3 Ni + 2 N b - 3, 

where N t is the total number of triangles and N e is the total number of edges. The impor- 
tance of these formulas (not so much what the values in the formulas are, but more the fact 
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that some fixed formula holds) will show up in the next section. If we let Af< represent the 
number of points joining to pi then it is easy to see that 

N 

]T Mi = 2 Ne 

t= 1 


and so we have that the “average valence” of a point is given by 


M = 


TaLi M i _ 6 0 ^ + 3 
Ni + N b N 


which is approximately 6. For a sphere (or any domain homeomorphic to a sphere) we 
have no boundary points and so N = Ni and the analogous formulas are 

g 

j\T t = 2[N — 1), N e = Z(N-l), M — -(N-l). 


Some Special Triangulations 

One of the simplest triangulations results from splitting the rectangles of a Cartesian grid. 
A Cartesian grid involves two monotonically increasing sequences, x t , i = 1, . . . , n and 
yj,j — 1, ... ,m. The grid points have coordinates (x, , j/j ) and these points mark out a 
cellular decomposition of the domain consisting of rectangles. See Figure 20. 10. Forming 
an edge with one of the diagonals of these rectangular cells leads to a triangulation of the 
domain. In Figure 20. 1 1 is shown a triangulation where a consistent choice for the diagonal 
is made. In Figure 20.12 is shown a triangulation with mixed choices for the diagonals. 
In some applications where dependent ordinate values are known, it is possible to base 
the choice of the diagonal upon some criteria such as minimum jump in normal vector 
(see Section 20.2.4) or whether or not the diagonal vertices are separated or connected 
based upon the hyperbolic contours at the mean value (see the asymptotic decider criteria 
discussed in [186]). In general for this type of triangulation which results from a Cartesian 
grid, it is not necessary to maintain the triangular grid structure (see Figure 20.8) as this 
information can be directly inferred from the natural labeling of p,j = (xi , y 3 ) . Only the 
information which indicates which diagonal is selected needs to be made available. 

We now want to discuss some special triangulations which result from curvilinear grids. 
A curvilinear grid is specified with two “geometry arrays” (x,j , Vij ) , i = 1 , - • • . M\ j = 
1 , . . . , jV. A cell Cij consists of the quadrilateral with the boundary delineated by (x,y , ) 

to {xi+n, yi+ij) to {x ij + u Vij+i) back to (xij, yij ). It is assumed that these four points 
form a simple (nonintersecting) polygon so that the quadrilateral is actually well-defined. 
This condition obviously puts some geometric constraints on the values of the geometry 
arrays that specify a curvilinear grid. 

An example of a curvilinear grid is shown in Figure 20.13. In this case the cell C 73 
degenerates to a triangle because ( As 3 > ^ 83 ) and ( , 184) are the same point and the cell 
Cs3 degenerates to an edge because, in addition, (A 93 , 193) and ( A94 , V’94) are the same 
point. The cells C 33 , C 43 , C 53 , C 63 , and C 73 have been removed from the domain creating 
the hole in the interior. 

The domain (the union of all of its cells) can be triangulated simply by triangulating 
each of the cells, by choosing a diagonal to an edge of the triangulation. An example related 
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Figure 20.13: An example of a curvilinear grid. 
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Figure 20.14: Triangulation resulting from curvilinear grid. 


to the grid of Figure 20.13 is shown in Figure 20.14. Here we have modified the grid by 
moving the point (X 72 , ¥ 72 ) a little. This serves to point out that if the cell is not convex, 
then there may be only one choice for the diagonal. 

We now discuss some special triangulations obtained by subdividing an existing trian- 
gulation. We briefly mention a couple of possibilities. The first is based upon inserting an 
additional point into the interior of an existing triangle and thereby forming three new tri- 
angles. This is illustrated in Figure 20.15. This particular type of subdivision is sometimes 
referred to as the Clough-Tocher split because of its association with a very well known 
finite element shape function defined over a triangular domain. 

Another way to subdivide an existing triangulation is to insert a new point on an existing 
edge and split the two triangles (unless the edge is on the boundary) which share this edge. 
If all edges are split simultaneously we obtain yet another triangulation where each previous 
triangle is replaced by four new ones. Two different ways for forming triangles from these 
points is shown in Figure 20. 16 and Figure 20.17, respectively. These types of subdivision 
are particularly interesting due to the nested properties of function spaces which are defined 
in a piecewise manner over the embedded subdivisions. This can lead to wavelets and their 
related multiresolution analysis. For the efficient application of these triangulations, it is 
important to have a method of labeling the triangles which allows an efficient algorithm 
for finding the labels of all neighbors of a triangle. The labeling scheme illustrated in 
Figure 20.17 has these properties. We call it the divide and flip scheme and have found it 
to be very useful for implementations. It is related to the spherical quadtrees discussed by 
Fekete [85]. The first step of the subdivision applied to the triangulation of Figure 20.8 is 
shown in Figure 20.18. 



Figure 20.16: Nested subdivision triangulation. 
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Figure 20.17: The divide and flip labeling scheme for a nested subdivision triangulation. 



Figure 20.18: A triangulation obtained by splitting each edge of an existing triangulation 
and forming triangles as indicated in Figure 20.17. 
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skinny triangle with 
a very small angle 


skinny triangle with 
a very large angle 


Figure 20.19: Examples of poorly shaped triangles. 


20.2.2 Optimal Triangulations 

Types and Characterizations 

There are many possible triangulations of a given, polygon-bounded domain D. For some 
applications (but not all) it is desirable to avoid poorly shaped triangles. See Figure 20. 19. 
These are triangles with very large angles or ones with very small angles. This gives rise to 
two types of optimal triangulations which have been discussed quite widely: the MaxMin 
and MinMax. Both of these optimal triangulations have a similar method of characteri- 
zation. Associated with each triangulation there is a vector with N t entries representing 
either the largest or smallest angle of each triangle. The entries of each vector are ordered 
and then a lexicographic ordering of the vectors is used to impose an ordering on the set of 
all triangulations. In the case of the MinMax criterion, A, is the largest angle of a triangle 
and the entries of each vector, A t , are ordered so that 

At — (Ai . • • • » Apff ) , Ai > Aj , i < j. 

The smallest of these vectors, based on their lexicographic ordering, associates with the 
optimal triangulation. In the case of the MaxMin criteria, a,- is the smallest angle and the 
entries of each vector are ordered the other way so that 

&t — (til i &2i ■ • • > a n , ) i a x ^ a j i * ^ j ■ 

The largest of these vectors represents the optimal triangulation in the MaxMin sense. In 
Figure 20.20, six data points are shown which have a total of ten possible triangulations, 
which are shown in Figure 20.21. The associated vectors for the MinMax criterion are 

A Ta = (2.84,2.36,1.99,1.77,1.57) 

= (2.98,2.84,1.99,1.91,1.57) 

Ar, = (2.98,2.42,1.91,1.88,1.57) 

A Ti = (2.84,2.36,2.32,1.99,1.40) 

A T4 = (2.42,2.36,1.88,1.77,1.57) 

A Ti = (2.98,2.42,1.95,1.91,1.27) 

A ra = (2.42,2.36,2.32,1.88,1.40) 

A T7 = (2.42, 2.36, 2.32, 1.50, 1-50) 

A Ta = (2.42, 2.36, 1.95, 1.74, 1.50) 

A T t = (2.42,2.36,1.95,1.77,1.27) 
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which we rearrange into decreasing order to obtain 

A Tl = (2.98,2.84,1.99,1.91,1.57) 

A t , = (2.98, 2.42, 1.95, 1.91, 1.27) 

A T2 = (2.98, 2.42, 1 .91, 1.88, 1-57) 

A t 3 = (2.84, 2.36, 2.32, 1.99, 1.40) 

A t „ = (2.84, 2.36, 1.99, 1.77, 1.57) 

A Te = (2.42,2.36,2.32,1.88,1.40) 

A T7 = (2.42, 2.36, 2.32, 1.50, 1 50) 

A r: , = (2.42, 2.36, 1.95, 1.77, 1.27) 

A r , = (2.42, 2.36, 1.95, 1.74, 1.50) 

A u = (2.42, 2.36, 1.88, 1.77, 1.57) 

which implies the following ordering 

U < T& < r<) < T7 < r$ < To < r 3 < t 2 < T S < ~i 

and so r 4 is the optimal triangulation in the MinMax sense. On the other hand, the associ- 
ated vectors for MaxMin criteria sorted in increasing order are 

a T , = (0.02,0.04,0.35,0.46,0.50) 

a T2 = (0.02, 0.11, 0.42, 0.46, 0.50) 

a T 5 = (0.02, 0. 1 1, 0.50, 0.58. 0.88) 

a Ti = (0.04, 0.14, 0.35, 0.37, 0.66) 

« r „ = (0.04,0.14,0.35,0.46,0.62) 

a T , = (0.11,0.14,0.37,0.42,0.66) 

a Tr = (0.11,0.14, 0.37,0.46,0.70) 
a T4 = (0.1 1, 0.14, 0.42. 0.46, 0.62) 

a Ts = (0.11, 0.14,0.57, 0.58,0.70) 

a T9 = (0.11, 0.14, 0.58, 0.62, 0.88) 

which results in the following ordering 


n < r 2 <t 5 < t 3 < r 0 < t 6 < r 7 < t 4 < r s < r 9 

and so r 9 is the optimal triangulation in the case of the MaxMin criterion. 

In the case where D is the convex hull of the points of P, there is an important rela- 
tionship between the MaxMin triangulation and the Dirichlet tessellation. The Dirichlet 
tessellation is a partition of the plane into regions R, , i = 1 , . . . , N called Thiessen re- 
gions. The Thiessen region R* consists of all points in the plane whose closest point among 
Pi , i = 1 , . . . , n is pk ■ A Dirichlet tessellation is usually illustrated by drawing the bound- 
aries of the Thiessen regions. The collection of these edges is sometimes referred to as the 
Voronoi diagram (see [252]). An example is shown in the left image of Figure 20.24. In 
the right image of Figure 20.24 is shown the MaxMin triangulation, which is also called 
the Delaunay triangulation. It is dual to the Dirichlet tessellation in that the edges of this 
optimal triangulation join vertices which share a common Thiessen region boundary. We 
have included the great circles in the left image of this figure so as to point out another 
important property of the Dirichlet tessellation and its companion Delaunay triangulation. 
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By definition, the edges of the Thiessen regions meet at triads (possibly more than three 
edges meet in some special, neutral/cyclic cases) which are equally distant to three points. 
These three points will form a triangle of the optimal triangulation and the great circle will 
not contain any other data points. See Figure 20.25. 

We can be a little more formal about these properties if we introduce some notation. 
Recall that I t = {(*(m) J< /(m), Ar(m)), m = 1,... , N t } so that the three data points 
Pi{m) » Pj(m) i Pk(m) will be the vertices of a triangle of the triangulation. We assume that 
the neighbor information of the triangular grid is given by three arrays ij(m ) , jk(m), and 
ki{m),m = 1, . . . , N t . Let V m be the point which is equidistant from and 

Pk(m) and C m = {p : ||p- V m \\ < \\V m - p a (m)\ \,a — h j or k} be the circumcircle (disk) 
for this triangle which has V m as its center. The Delaunay triangulation is characterized 
by the fact that C m does not contain any other data points p^ i = 1, . . . , N other than 
Pi(m) ? Pj ( m ) j and Pk(m) • The points V m are the vertices of the Voronoi diagram. In order 
to draw the Voronoi diagram we simply start with some V m and draw the edges to the three 
points that are joined to it, namely, Kj(m)i Vjk{m)} Vki{m)> If any one of ij(m) 1 jk(m) } 
or ki(m) is zero (say ij(m ), indicating the edge joining p l(m) and p j(m) is on the boundary 
of the convex hull), then we draw the ray emanating from V m in the direction perpendicular 
to the appropriate edge (which isp i{m )p j(m) if ij(m) = 0, p j{m) p k{m) if jk(m) = 0, and 
Pk{m)Pi{m) ^ ki(m) — 0). If we go through the list of triangles and draw three edges for 
each V m we will actually be drawing each edge (not each ray) twice. We can avoid this 
duplication by testing (for example) whether or not m > ij(m), m > jk(m ), m > ki(m) 
before we draw the corresponding edge. 
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Figure 20.23: Spherical triangulation and tessellation. 


Because of this relationship between the Dirichlet tessellation and the optimal MaxMin 
tri angulation, we can extend the idea of MaxMin or Delaunay triangulation to any domain 
where we can compute the distance between two points. The sphere provides an interesting 
and useful example. Here the distance between two points p and q is easily computed as 
cos - 1 (p * q) so the Dirichlet tessellation is also easy to compute. An example is shown in 
the right image of Figure 20.23. The left image depicts the triangulation which is dual to 
this tessellation. 

There have been many other criteria for characterizing optimal triangulations that have 
been studied and discussed in the literature. Some turn out to be equivalent to those we 
have mentioned here and some only appear to be similar, so one needs to be rather careful. 
Even though the terminology can be similar, the criterion of minimizing the maximum 
angle is not the same as the MinMax criterion we have described here. It is easily the 
case that the two quite different triangulations with different vectors A T (as defined above) 
could have the same maximum angle and could both be a triangulation which minimizes the 
maximum angle. The example of Figure 20.20 has this property. Each of the triangulations 
r 6 , ty, r 9 , rs, and r 4 have a maximum angle of 2.42, which turns out to be a minimum, 
so any one of these triangulations would satisfy the criterion of minimizing the maximum 
angle, while only r 4 satisfies the MinMax criterion described here. Overall, the topic of 
optimal triangulations can be rather technical, and one has to be careful when comparing 
results found in the literature. 
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Figure 20.24: The Dirichlet tessellation and its dual triangulation. 

Algorithms for Delaunay Triangulations 

In this section we discuss some ideas and techniques leading to algorithms for computing 
the Delaunay triangulation of a set of points in the plane. In general, this is a very rich 
and full area of research and here we can only provide a glimpse. The literature is very 
abundant with both practical and theoretical papers on this subject. There is not a single 
“best” algorithm. The choice depends upon the particular application and the tools and 
resources available. It is a good strategy to be armed with a collection of ideas, tools, 
and techniques so that an effective algorithm can be custom-designed for the application at 
hand. Our approach for the material for this section is based upon a discussion of the ideas 
behind a few selected algorithms. Our selection is based upon potential usefulness of the 
ideas and on which would be representative. In addition, we are particularly interested in 
those ideas which extend most easily to three dimensions. But, just for the sake of interest, 
we have included the description of one 2D algorithm which does not extend at all to 3D! 

The Swapping Algorithm of Lawson [139]: The basic operation of this algorithm 

consists of swapping the diagonal of a convex quadrilateral. Lawson [138] showed that any 
triangulation of the convex hull can be obtained from any other triangulation by a sequence 
of these operations. (Later this property was established for nonconvex domains by Dyn 
and Goren [66].) Furthermore, Lawson proved that if the choice of the diagonal is made on 
the basis of the MaxMin criterion for the quadrilateral only, eventually the global optimal 
triangulation will be obtained. In other words, for this criterion, a local optimum is a global 
optimum. A typical implementation of this type of algorithm would insert new points (say, 
in sorted x-order) in the interior of an existing triangulation or connect to all points on the 
boundary which are visible from the new point. This new triangulation is then optimized 
by testing and possibly swapping the diagonals of convex quadrilaterals. It is interesting to 
note that this type of algorithm will not necessarily produce the MinMax because for this 
criterion, a local extreme is not necessarily a global optimum. The example of Figure 20.20 
illustrates this. Based upon the MinMax criterion, r 4 is optimal and r$ is a local minimum. 
Locally optimal swaps of diagonals from r 8 would never lead to r 4 . The algorithm could 
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Figure 20.25: Notation and terminology for Delaunay triangulation and Dirichlet tessella- 
tion. 
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easily get trapped in a local extreme at r 8 . The ideas of simulated annealing can be used to 
develop algorithms which can escape from these local extrema. See Schumaker [225], for 
example. 

The Algorithm of Green and Sibson [107]: This algorithm depends heavily upon a 
particular data structure used to store the Delaunay triangulation (or Dirichlet tessellation). 
For each object (a Dirichlet tile or window boundary constraint) is recorded in a “contiguity 
list” consisting of all objects with which it is contiguous. This data structure is very similar 
to the contiguity list structure we described in Figure 20.9 but it also includes some window 
boundary constraints. New points are inserted sequentially. We quote directly from [ 1 07] 
to how this is done. 

The contiguity list for the new point is then built up in reverse (that is, clock- 
wise) order and subsequently standardised. We begin by finding where the 
perpendicular bisector of the line joining the new point to its nearest neigh- 
bour meets the edge of the nearest neighbour’s tile, clockwise round the new 
point. Identifying the edge where this happens gives the next object contigu- 
ous with the new point and this is in fact the first to go onto its contiguity list. 

The new perpendicular bisector is then constructed and its incidence on the 
edge of this new tile is examined to obtain the subsequent contiguous object; 
successive objects are added to the contiguity list in this way until the list is 
completed by the addition of the nearest neighbour. Whilst this being done 
old contiguity lists are being modified: the new point is inserted in each and 
any contiguities strictly between the entry and exit points of the perpendicular 
bisector are deleted, the anticlockwise-cyclic arrangement of the lists making 
both this and the determination (sic) of the exit very easy. 

This insertion algorithm requires the computation of the nearest existing data point to the 
data point that is to be inserted. The authors discuss an algorithm which takes advantage of 
the tessellation computed so far. In the authors’ words: “Simply start at an arbitrary point 
and “walk” from neighbour to neighbour, always approaching the new point, until the point 
nearest to it is found.” 

The Algorithm of Bowyer [21]: Bowyer described an algorithm for inserting a new 

point (lying in the convex hull) into an existing Delaunay triangulation. An example given 
by Bowyer and which we include in Figure 20.27 serves to define this data structure. In 
the terminology of Bowyer, the forming points for a vertex are simply the vertices of the 
triangle which has this particular vertex as the center of its circumcircle. Since each triangle 
gives rise to a vertex, giving a list of indices of the forming points for each vertex (as 
Bowyer does) is equivalent to giving a list of indices of the data points which comprise 
each triangle of trianguiation. Except for a change in ordering, the neighboring vertices are 
exactly the same as the indices of the triangle neighbors as given in the triangular grid data 
structure of Figure 20.26. 

In order to insert a new point ( Q in Figure 20.27) within the current convex hull of the 
data points, Bowyer [21] gives the following algorithm: 



20.2 Triangulations 


441 



Figure 20.26: An aid to the Green and Sibson algorithm. 

1 . Identify a vertex currently in the structure that will be deleted by the new point (say 
V 4 ). Such a vertex is any that is nearer to the new point than to its forming points. 

2. Perform a tree search through the vertex structure starting at the deleted vertex look- 
ing for others that will be deleted. In this case the list will be: {V 4 ,V 3 ,V 5 }. 

3. The points contiguous to Q are all the points forming the deleted vertices. 
{Pa, ft. A, ft, *r}. 

4. An old contiguity between a pair of those points will be removed (ft — P 4 say) if all 
of its vertices { V^, ft) are in the list of deleted vertices. 

5. In this case the new point has five new vertices associated with it: 
[Wx , W 2 , W 3 , W 4 , W$}. Compute their forming points and neighbouring vertices. 
The forming points for each will be the point Q and two of the points contiguous to 
Q. Each line in the tessellation has two points around it (the line V 3 - V 2 , for example, 
is formed by P 3 and ft). The forming points of the new vertices and their neigh- 
bouring vertices may be found by considering vertices pointed to by members of the 
deleted vertex list that are not themselves deleted, and finding the rings of points 
around them. Thus If 5 points outwards to V 2 from Q and is formed by { P 3 ■ ft) Q } • 

6. The final step is to copy some of the new vertices, overwriting the entries of those 
deleted to save space. 

The Algorithm of Watson [254]: This algorithm relies on the property of a Delaunay 

triangulation that a triple of data point indices (i,j, k) will be in I t provided the circum- 
circle of pi,Pj, and p k contains no other data points. As with the other algorithms, this 
algorithm is based upon inserting a new point into an already existing Delaunay triangu- 
lation. The general philosophy of Watson’s approach is described by the following two 

steps: 
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1 . Find all triangles whose circumcircle contains the point to be inserted. 

2. For each of these triangles, form three new triangles from the point to be inserted 
and the three edges of this triangle and test to see if any of these three new triangles 
contain any other data points. If not, then add this new triangle to the triangulation. 

More details for this general approach are given in the flow diagram of Figure 20.28, 
which is based upon the flow diagram of [254], 

Watson [254] describes a number of features and details to make the basic algorithm ef- 
ficient and eventually discusses a particular implementation which he says has an expected 
running time which is observed to increase not more than N 3/2 . 

The Embedding/Lifting Approach: Algorithms of this type are based upon a very 

interesting relationship that exists between the three-dimensional convex hull of the lifted 
points (xi,yi,xj + yf) and the Delaunay triangulation. Faces on the convex hull are des- 
ignated as being either in the upper or lower part. The lower part consists of faces which 
are supported by a plane that separates the point set from (0, 0, -oo). The Delaunay tri- 
angulation is obtained directly from the projection onto the £ (/-plane of the lower part of 
the convex hull. See [27] and [68]. An algorithm for computing the convex hull which 
is based on an initial sort followed by a recursive divide-and-conquer approach has been 
described by Preparata and Hong [202]. This algorithm is also covered in [68] and [203]. 
Theoretically the algorithm is optimal time O (n log(n)), but Day [49] reports that empir- 
ical data implies a worst-case complexity of 0(N 2 ). Day covers many of the details and 
special-case issues of practical interest for implementation which are often brushed over in 
more theoretical papers. 

Divide-and-Conquer Algorithms: The general structure of this type of algorithm is 
to divide the data set into subsets A and B, solve the problem for A and solve the problem 
for B and merge the results into a solutionfor A U B. See Figure 20.29. Divide-and-conquer 
algorithms can lead to theoretically optimal algorithms, but often fail to be competitive in 
practical usage. The merging portion is often the most troublesome in trying to maintain 
bounds on the running times and complexity of the algorithm. 


20.2.3 Visibility Sorting of Triangulations 

This is an example of an area that is interesting in 3D but not in 2D. It is possible to make a 
definition of a visibility sort for a triangulation which is completely analogous to that of a 
tetrahedrization, but there does not appear to be any application or use for such a property. 
We defer further discussion on visibility sorting to Section 20.3.3. 

20.2.4 Data-Dependent Triangulations 

The topic of data-dependent triangulations arises within the context of determining a mod- 
eling function F(x, y) for the data (F, ; £,, y,), i = 1, . . . , N. A relatively simple approach 
to defining a modeling function is to first form a triangulation of the convex hull of the 
independent data (x< , yj), i = 1 , . . . , N and then define F to be piecewise linear over this 



444 


Tools for Triangulations and Tetrahedrizations 



Figure 20.28: Flow diagram for Watson’s algorithm. 
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Figure 20.29: Divide-and-conquer algorithms. 




triangulation. This will yield a C° (continuous) function which interpolates the data; that 
is, F(xi,yi) = Fi,i = 1, . . . , N. We denote this function by F T (x, y). Any triangula- 
tion of the independent data (x it »), * = 1, • • • , N will suffice for this approach. While 
we are well aware of the many desirable properties of the Delaunay triangulation, it might 
very well be the case that some other triangulation whose choice would depend upon the 

values F{J = 1 ,N would lead to some desirable properties for the modeling function 

F. This is the basic idea of data-dependent triangulation. Of course, there are potentially 
many ways to accomplish this, but we choose for this discussion here to briefly describe 
the criteria called “nearly C l ” as proposed in [67]. An ordering is imposed on the collec- 
tion of all possible triangulations of the convex hull in the following manner. First a local 
cost function (see Figure 20.30) for each edge e,- = 1, . . . , N{ e = N e — Nf, is defined and 
denoted by S{ F T , e { ). (We will shortly describe the four examples of local cost functions 
covered in [67].) If T and T f are two triangulations, then 

T <T' 

provided the vector 

(s(Fr,ei), s(Fr,e 2 ), * • • ,*(*r,<ov, e )) 
is lexicographically less than or equal to 

{s(F T 'yei),s(FT^e 2 ), . . • , s{Ft', eN te ))■ 

It is assumed that the components of these vectors are arranged in nonincreasing order. The 
goal is then to find the optimal data-dependent triangulation which is defined by having 
the smallest associated vector under this lexicographical ordering. Since there are only a 
finite (albeit possibly very large) number of possible triangulations, we know that a global 
minimum exists even though it may not be unique and it may not be so easy to compute. 
The algorithm used in [67] is similar to the swapping algorithm of Lawson (which we have 
described above in Section 20.2.2) in that an initial triangulation is obtained and then an 
internal edge of a convex quadrilateral is considered. If V < T, where V is the same 
triangulation as T except the diagonal of the convex quadrilateral has been switched, then 
this switch is made and other edges are considered for potential swapping. Since each swap 
moves strictly lower in the lexicographic ordering, we are guaranteed that this algorithm 
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Figure 20.30: Notation for local cost function definitions. 


will eventually converge after a finite number of steps. This means that swapping any 
edge would not move to a smaller triangulation. This limit triangulation may not be the 
global minimum; it is only guaranteed to be a local minimum and steps to find the global 
minimum must do more than swap diagonals which improve (with respect to the ordering) 
the triangulation. 

We now describe the four local edge cost functions used in [67]. Let P\ = aix+b\y-\-ci 
and Po — a 2 x 4- b 2 y 4- co be the two planes defined over the two triangles of a convex 
quadrilateral. 

i) The angle between normals: The local cost function is taken as the acute angle be- 
tween N\ and iV 2 , which are the respective normals for P\ and P 2 . 

s(Fr , e) = cos _1 (/l) 


where 


q _ a i a 2 4- b x b 2 4* 1 

V^( a i TW^T) [a\ 4- 4- 1) 

ii) The jump in normal derivative: This cost function is the difference between the 
derivative of Pi and P 2 - This derivative is taken in the direction perpendicular to 
the edge dividing the two triangles. 

s(F T ,e) = [n x (a x - a 2 ) + n y {bi - i 2 )] 
where (n T , n y ) is a unit vector perpendicular to the edge e. 

iii) The deviations from linear polynomials: The cost function measures the error be- 
tween Pi and P 2 , evaluated at the other point of the quadrilateral. 


s(F T .e) - \fFF~yi) - ^.) 2 + (Po (x k ,y k ) - Ff 
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Figure 20.31: Examples of data-dependent triangulations. 


iv) The distance from planes: This cost function measures the distance between the 
planes P\ and P 2 and the corresponding vertex of the quadrilateral. 


/r , x [{Piixi/yi) - Fi ) 2 {P2(xk,yk) - F k f 

s(Ft ' e) = V a? + 6?+l"' + ^T+fcT+l 

Some typical results are given in [67] which confirm the expectation that using the 
optimal data-dependent triangulation improves the overall fitting properties of F T over that 
of the Delaunay triangulation, which, by the way, is used as the initial triangulation for the 
swapping algorithm. It is observed that long, thin triangles tend to appear where the data 
seems to indicate a function that is increasing (or decreasing) relatively rapidly in a certain 
direction. The use of the data-dependent triangulation generally gives an overall reduction 
in errors when certain test functions are used to generate the data. 

As we have mentioned, the local swapping algorithm used in [67] can only find a local 
minimum. In order to move more closely to the globally optimal data-dependent triangula- 
tion, Schumaker [225] and Quak and Schumaker [204,205,206] have involved the tools of 
simulated annealing. More details on this are contained in Section 20.3.6 on data-dependent 
tetrahedrizations. We include here the results of one example described by Schumaker. The 
data consists of 


(Fij'.Xi, yj); Xi,yj 


0.0, 0.2, 0.4, 0.6, 0.8, 1.0; 


where 


F{x,y ) = (y- z 2 )+- 

Three triangulations are shown in Figure 20.3 1 . The first is the Delaunay triangulation of 
the independent data. The next is the triangulation which results from the local swapping 
algorithm of [67] using the local cost function of “angle between normals.” The last is the 
triangulation after simulated annealing has been applied. The associated vectors for each of 
these triangulations is given in Figure 20.33 and the piecewise linear approximations over 
these triangulations are shown in Figure 20.32. 
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Figure 20.32: The graphs of Schumaker’s example. See [225]. 
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Angles between normals for Delaunay triangulation: 
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Figure 20.33: Angles for the data- dependent triangulation. 
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Figure 20.34: Two different units used to measure the same data lead to two different 
Delaunay triangulations. 


20.2.5 Affine Invariant Triangulations 

The desirable properties of the Delaunay triangulation have been previously discussed. Un- 
fortunately, this optimal triangulation is not invariant under affine transformations, and this 
means that methods for analyzing and visualizing data that use this particular triangulation 
can be affected by the choice of units used to measure the data. This could be considered an 
undesirable property. In this section we describe a relatively new method for characteriz- 
ing and computing an optimal triangulation which is invariant under affine transformations. 
Before we proceed with the discussion of these techniques, we wish to motivate further the 
desirability of affine invariance. 

As we have mentioned earlier, one of the main purposes for triangulations and tetra- 
hedrizations is their use in defining functions in a piecewise manner over the domain of a 
data set. It would be undesirable if the happenstance of the choice of units used to measure 
the data were to affect the definition of a data modeling function. But this does happen 
with the Delaunay triangulation. The example of Figure 20.34 points this out. This data 
represents the independent data; the dependent data is not given as it is not important in 
this context. The data is the same in both the left and right graphs of Figure 20.34; the only 
difference is that in the left graph we have used years and £ (pounds, British monetary 
unit, approximately equal and assumed here to be exactly equal to two US dollars), and 
in the right graph we have used months and dollars. If we use the units of years and £ 
then we can see that the three vertices (lyr, \£), (0.5yr, 3£) y (2yr, 2£) will mark out a 
triangle to be included in the list of triangles for the Delaunay tri angulation. But on the 
other hand, if we use months and $ we can see that the circumcircle defined by these same 
three vertices (12mon, 2$), (6mon, 6$), (24mon, 4$) contains the data point (4mon, 4$). 
Therefore, these three vertices will not comprise a triangle of the Delaunay triangulation if 
these units are used. This simple example points out the possible effects of the choice of 
the units of measurement. The choice of the units of measurement is the same as a change 
in scale, x ax and y <— by. Uniform scale changes of the type x ar, y ay will 
not affect the Delaunay triangulation. 

We now discuss how to avoid this problem. It would be possible to simply normalize 
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all data ranges to one unit by scaling by the range. But this approach would mean that 
rotations of the data could have an effect on the Delaunay triangulation, meaning the final 
data model would be affected by rotations of the data. In other words, the placement and 
alignment of the axes for the measurement of the data would have an effect on the data 
modeling function and subsequently on our analysis of the data, and this we would like to 
have the opportunity to avoid. It would, in general, be useful to have a characterization 
(and subsequent algorithms) for an optimal triangulation which is not affected by affine 
transformation. An affine transformation is a map of the form 


(x,y) = A{x,y) + c 


where A is a 2 x 2 matrix and c is a two-dimensional point. Affine transformations in- 
clude not only scale changes and rotations, but also translations, reflections, and shearing 
transformations. The approach to such an optimal triangulation covered here is through 
the duality that exists between the conventional Delaunay triangulation and the Dirichlet 
tessellation. As we described previously, the characterization of the Delaunay triangulation 
(as a MaxMin triangulation), it was heavily dependent upon angles, and angles are affected 
by scaling transformations; so it should be no surprise that the Delaunay triangulation is 
also affected by scaling transformations. But the definition of the Dirichlet tessellation 
uses only distance and we know that the Delaunday triangulation is dual to (a direct result 
of) the Dirichlet tessellation. The approach here is to use a method of measuring distance 
which is invariant under affine transformations. The Dirichlet tessellation based upon this 
new method of measuring distance will have a dual which will serve as our op timal trian- 
gulation. Rather than use the standard Euclidean norm | j (ar , j/)||" = \fx - 4- j/ 2 we propose 
the use of the following norm 


||(z,2/)llv = ( x <y) 


where 
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We have used the subscript of V on the norm to explicitly indicate that this method of 
measuring distance is dependent upon the data set. Change the data set and you change how 
you measure distance, but the distance between any two data points will remain constant. 
This norm and its use within the context of scattered data modeling was first described in 
[181], This norm has the property that it is invariant under affine transformations. More 
precisely. 


\\P-Q\W =\\T{P)-T{Q)\\ t[V) 

for any two points P = (x, y) and Q = (u, v) and any affine transformation 
Here, T(V) (used as a subscript in Equation (20.2)) is the transformed data 


( 20 . 2 ) 


T[V) = ( T[ Xl ^ 

'2/1- Py 


*2 “ fix 
Vi ~ Vy 


T 


•i\V ~ fi x \ \ 

yy ~ fly ) ) 


Figure 20.35 illustrates the properties of this new method of measuring distance. Each 
of the data sets shown in this figure are affine images of each other. Starting in the upper left 
and moving in a clockwise direction, the transformations are: counterclockwise rotation of 
44 degrees; a scaling in x by a factor of 2; a scaling in y by a factor of 0.4. The four 
ellipses in each figure represent points which are 1/4, 1/2, 3/4, and 1 unit(s) from their 
center point as measured with the affine invariant norm. In Figure 20.36 we show the 
Dirichlet tessellation of these four affinely related data sets, and in Figure 20.37 we show 
the corresponding dual triangulation. As one can see, the triangulation is unchanged by 
these transformations. 

As a comparison, we have also included the Delaunay triangulation based upon the 
standard Euclidean norm in Figure 20.38. And as we indicated earlier, we can see that 
triangulation results are affected by the transformations. Not all triangles are changed, but 
some are. 


And now we suggest some practical information on how to incorporate this feature into 
an algorithm for computing triangulations. If you already have a procedure for computing 
an optimal triangulation, then it is possible to modify it slightly to achieve the results we 
have described in this section. Say, for example, that the procedure is based upon Law- 
son’s algorithm and that there is a subprocedure which decides whether or not to switch 
the diagonal of a quadrilateral formed from two triangles. It might be that this procedure 
is based solely on Euclidean distance. That is, the center and radius of the circumcircle of 
three points are determined and the distance to the center from the fourth point is computed 
so as to make this decision. In order to modify this subprocedure, we need only to replace 
the use of the Euclidean norm with the affine invariant norm described here. The equations 
for computing circumscribing circles (ellipses) for a quadratic norm in general are given 
in [182]. If, on the other hand, the procedure you are already using is known to be rota- 
tion invariant, then there is an even easier way to affect the results of the affine invariant 
triangulation. This is based upon the factorization of the matrix which defines the affine 
invariant norm. We denote this matrix by A{V) so that we have 


ll(*.y)llv = ^ * ) 
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The matrix A{V) can be factored (Cholesky) into 

^)=(!;i £)('» t)=ww- 

Here the notation L(V)‘ denotes the transpose of L(V). 

Using this factorization, we have that 

||(x,y)||t- = (x,y)L(V)L(Vy (*) = IK*. y)W )« 2 

which means measuring distances with the affine invariant norm is the same as measuring 
distance in the standard Euclidean but with the points transformed by multiplying by L(V). 
This means that we can achieve the result of the optimal affine invariant triangulation by 
computing the standard Delaunay triangulation on the transformed data 

(X i ,Y i ) = (x il y i )L(V) 


In summary, we need only to compute 
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to apply any rotation invariant triangulation algorithm to the transformed data 
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20.2.6 Interpolation in Triangles 

We now take up the topic of interpolating into (or over) a single triangular domain. The 
interpolants we describe here form the basic building blocks for constructing the global in- 
terpolants which have piecewise definitions over the individual triangles of a triangulation. 
The domain here is a single triangle, T = Tijk with vertices Vj,Vj, and Vk , and the data 
consists of values given on the boundary of the triangular domain. We need to differentiate 
between two types of boundary data. If the data consists of function and certain derivative 
values specified only at the vertices (or possibly other points such as midpoints), then we 
call this discrete data. If, on the other hand, the data is provided on the entire boundary 
of the triangle, we refer to this type of data as transfinite data. The importance of an in- 
terpolant which will match transfinite data is that it serves as a prototype for developing 
a large variety of discrete interpolants. This is accomplished through the process of dis- 
cretization, where the data required for a transfinite interpolant is provided by means of 
using some interpolation scheme only on the boundary, discrete data. For example, given 
only data values at the vertices, we can use linear interpolation along an edge to produce 



458 


Tools for Triangulations and Tetrahedrizations 


the transfinite data required by the transfinite interpolant, or if we also have data values at 
the midpoints, we could use quadratic interpolation. 

There is a second concept which is rather important for interpolants defined over trian- 
gles and this has to do with the degree of continuity of the global interpolant. Often, we 
require that the global interpolant at least be continuous. We call such an interpolant a C° 
interpolant. If the global interpolant has continuous first-order derivatives, we say it is a 
C 1 interpolant. A C° interpolant for a single triangle is one which interpolates to boundary 
data consisting of only position values, either at the vertices only (and possibly points along 
the edges) or on the entire boundary. A C 1 interpolant for a single triangle is one which will 
interpolate to first-order derivative data specified on the boundary. But this must be done in 
a manner so as to guarantee C 1 continuity across the boundary edges. So, for example, if 
the cross-boundary derivative varies quadratically along an edge, then the data on this edge 
must be sufficient to uniquely determine this derivative, so that on an adjoining triangle we 
will have exactly the same cross-boundary derivative. For this reason, it is common for C 1 
interpolants to have linearly varying cross-boundary derivatives which are determined by 
their values at the two endpoint vertices. 

Combining the two concepts of discrete and transfinite data and C° and C 1 data leads 
to four types of triangular interpolants. This general area of interpolation in triangles is 
fairly rich and well developed, and we urge the really interested reader to follow the ci- 
tations into the literature (for example [177], [178], and [189]) after taking a look at the 
sampling we have chosen to include here. We first cover C°, discrete interpolants, then a 
sampling of three C°, transfinite interpolants. This is followed by the description of a C 1 , 
discrete interpolant. We have chosen to include a discretized version of the minimum norm 
triangular interpolant (see [178]). Another rather popular C\ discrete interpolant, is the 
Clough/Tocher interpolant often mentioned in conjunction with the finite element method. 
Much has been written about this interpolant in the past and so we do not include it here. 
This section is concluded with a description of a C\ transfinite interpolant, called the side- 
vertex interpolant [177]. It is one of the easiest to describe and the most versatile to use. It 
also generalizes rather nicely to a tetrahedral domain. 


C°, Discrete Interpolation in Triangles 

The lowest-degree polynomial, C°, discrete interpolant, is linear and unique. Given the 
data F(Vi) , F(Vj) , and F( V* ), the coefficients of the linear function 

F(x , y) = a + bx + cy 

which interpolates this data can be found by solving the linear system of equations 


a 4- bxi + cyi — F(VJ) 
a + bxj + cyj = F(Vj) 
a + bx k + cy k = F(V k ). 


Another path to this basic linear interpolant is via barycentric coordinates. Given a point 
V — ( x , y), barycentric coordinates, 6;, bj , and b k of this point relative to the triangle Tij k 
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are defined by the relationships 


x 

y 


= biV t + bjVj + b k V k 


= bi + bj -h b k ■ 


The linear interpolant now takes the form 


F(x, y) = F{ V) = biF(Vi) + bjF{Vj) + b k F(V k ). 


There are several alternative ways of defining or determining the barycentric coordi- 
nates of a point. For example, 


bi = 


A 

.4 




where A , , Aj , and A k represent the areas of the subtriangle shown in Figure 20.39 and A 
is the area of Tij k . Also, 
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X ~~~ Xj X | Xj 

y - t/j yi - yj 

X k “ Xj Xi Xj 

y k - yj yi - Vi 


Given the values at the three vertices and the three midpoints of a triangle, there is a 
unique quadratic which interpolates this data, 

Q(x,y) = F{Vi)bi(bi-bj-b k ) + F{Mj k )Abjb k 
+F(Vj)bj(bj -bi- b k ) + F(M ik )4bib k 
+F(V k )b k {b k - bi - bj) + F(M,j)4b,bj 

where M jk = ( V } + V k )/‘2, M, k — {Vi i + V k )/2 and M t] = (V x + Vj)/ 2. 

A common way to specify a cubic along an edge is to use the Hermite form which 
involves the first-order directional derivatives along the edges 

FUVi) = (x k - Xi)F x {Vi) + ( y k - Vi)F y {Vi) 

which are further illustrated in Figure 20.40. 
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Figure 20.41: The side- vertex interpolant notation. 


The six directional derivatives at the three vertices along with F(Vi), F{Vj) and F ( V k ) 
do not uniquely determine a cubic since the bivariate cubics are of dimension 10. The 
interpolant 

C(x,y) = F(Vi)b; + F^(Vi)b% + F; t (Vi)bjbj 
+F(V j )b 2 j + F!j[Vj)bjbi + F' k3 {V 3 )b]b k 
+F(V k )b k + F' ik (Vk)b k bi + Fj k (V k )b k bj 
+wbibjbk 


will match this function and derivative data for any value of w. This remaining degree of 
freedom represented by w can be absorbed by a variety of conditions. For example, it can 
additionally be required that the interpolant match some predescribed value at the centroid. 
Another common choice is 

w = 2[F(V5)+m) + F(l4)] 

+i[F^,(K) + Fji(Vi) + F-j(Vj) + Fij(Vj) + F'k(V k ) + F' k (V k )] 

which guarantees quadratic precision and is a result of discretization of a number of trans- 
finite interpolants (see [189]). Quadratic precision means that whenever the data comes 
from a bivariate quadratic function the interpolant will become this very same quadratic 
polynomial. 


C°, Transfinite Interpolation in Triangles 

In this section, we give only a sampling of three interpolants which will interpolate to arbi- 
trary function values on the boundary of a triangular domain, T = T ijk . More information 
on this general topic can be found in [189]. 
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Figure 20.42: The evaluation points (stencil) for side-side interpolant. 


The Side- Vertex interpolant: The side- vertex interpolant is built from three basic 
interpolants which are defined by linear interpolation along line segments joining a vertex 
and the opposing side. See Figure 20.41. In terms of barycentric coordinates, we have 


A[F] = biF(Vi) + (l- bi)F(S t ). 
MW = >> j F(V j ) + (l~b j )F(S J ), 
A k [F ] = + (1 - 6 fc )F(Sft) 


where Si 


Each of these interpolants will 


b J V l v * C. — C _ 

bj+bk “ 6 t +6 fc '° k “ b t +bj ■ 

interpolate to arbitrary function values on one edge of the triangular domain. In order to 
obtain an interpolant which matches arbitrary values on the entire boundary of T ijk , we 
form the Boolean sum of these three interpolants: 


A[F] = .4, © Aj © A k [F] = Ai [F] + Aj [F] + .4* [F] 

-Ai[Aj[F]] - Aj[A k [F]] - A k [Aj[F}} + .4,[.4 ; [4fc[F]]] 
= (1 - bi)F(Si) + (1 - bj)F(Sj) + (1 - b k )F(S k ) 
-biF(Vi) - bjF(Vj) — b k F(V k ) 


The Side-Side Interpolant: The side-side interpolant (Figure 20.42) is based upon the 
basic operation of linear interpolation along edges which are parallel to the edges of 
There are three of these interpolants: 


Pi[F] = 


Pj[F) = 


Pk[F] = 


b k F(bjVi + (1 - bj)V k + bjFjbjVi + (1 - 

b k + bj 

bjF^bjVj + (1 - bj)Vi + b k F(bjVj + (1 - bj)V k 

bi + b k 

biF(b k V k © (1 - b k )Vj + bjFjbkVk + (1 - b k )Vj 

bi + bj 



20.2 Triangulations 


463 



Unlike the basic interpolants of the side-vertex interpolant, these interpolants do not 
commute and so their triple Boolean sum is not well defined. However, it is possible 
to form the average of all double Boolean sums (each of which interpolates to the entire 
boundary) to arrive at the following affine invariant interpolant 


bkFjbjVi + (1 - bj)Vk) + bjF{biVi + (1 - bj)Vj) 

bk 4- bj 

bjFjbjVj + (1 - bj)Vi) + bkFjbjVj + (1 - 6j)Vfc) 

bi + bk 

bjFjbkVk + (1 - b k )Vi) + bjF^bkVk + (1 - b k )Vj) 

6 > + bj 

-biF{Vi) - bjF(Vj) - b k F(V k ). 
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The C* Interpolant: The third and final transftnite, C° interpolant we describe utilizes 
the stencil illustrated in Figure 20.43. 


C*[F](bi,bj,b k ) 


bib] 
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2b{ bj bk 


(b % 4- 2ij)(6j 4- 26, ) 
which can be written in the form 




C m [F](bi,bj,b k ) = b i F(V i ) + b j F(V j ) + b k F(V k ) 

+W k (f(Q k ) - (bi + b -fj F(Vi) - (bj + ^ F(Vj ) | 

+ ^ |f(Q,) - (ft, + ^ F(K) - (b k + ^ F(V'fc)} 

+Wi {m) - (bj + |) F ( Vj ) - (bk + F( 14)1 

where 

Wi = 4bjbk , l Vj = iM* , w k = 4b ' hj 

(2 bj + b t )(2bk 4- bi) (2 bi 4- bj)(2bk 4- bj) (26 t 4- 6^ ) (26j +6/;) 


Qi ~ ( bj + 1 ) Vj + ( bk + 1 ) 14 • 

Qj — ^ b{ 4- Vi + ^bk 4- 7 ^ 14 , 

Qc- = + y) V + + y) v j- 


In this form of C" we can see that it consists of linear interpolation plus a correction term. 
It can easily be verified that C* is precise for all quadratic functions. That is, if / is a 
quadratic, bivariate polynomial, then C*[/] — /. 
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C 1 , Discrete Interpolation in Triangles 

A commonly used 9-parameter, C 1 interpolant is 

C a [F](x, y) = {F{Vi) [&?(3-26<) + 6u>6i(6fca 0 + 6*0^] 

(* 

[6, r 6^ + wbi(3bitCiij + bj — 6^)] 

+Fji{ Vi) [bjbj + wbi(3bjOiik + b k - &,•)] } , 

where 

F' ki {Vi) = ( x k - Xi)F x {Vi) + {y k - yi)F y {V t ), 

Fji(Vi) = (xj - Xi)F x (Vi) + (y 3 - yi)F y (Vi), 


bibj b k 

bibj + b^k + bjbk 


I = {{i,j,k),(j,k,i),(k,i,j)}. 


and 

_ {Mf + H e *fcll~ ~ ll e u ll~ 

“ 2||c,-*|P 

We use \\eij\\ to denote the length of edge e ij; This 9-parameter, C 1 interpolant is a 
discretized version of a transfinite, C 1 , triangular interpolant, which is described in [178]. 
The derivatives which are in a direction perpendicular to an edge vary linearly along an 
edge. This guarantees that when two of these interpolants share a common edge the two 
surface patches will join with continuous first-order derivatives. It is possible to discretize 
the same transfinite interpolant and use an additional three parameters consisting of cross- 
boundary derivatives at the midpoints of the three edges. This leads to an interpolant that 
has all first-order derivatives vaxying quadratically along the edges. For a comparison of 
the Ca interpolant to the Clough/Tocher interpolant within the context of triangle-based 
scattered-data models, see Franke and Nielson [97]. 

C 1 , Transfinite Interpolation in Triangles 

In this section, we extend the problem of interpolating to transfinite data on the boundary 
to include also the requirement that the interpolant match user-specified transfinite deriva- 
tive data on the boundary. These types of interpolants can be used to construct surfaces 
over triangulated domains which are — that is, functions which have continuous first- 

order partial derivatives. One of the most versatile and easily described C , transfinite 
interpolants is the C 1 , side- vertex interpolant [177]. 

Earlier, we saw that the basic building blocks of the C°, side- vertex interpolant con- 
sisted of linear interpolation along lines joining a vertex and its opposing side. In order 
to extend these ideas to C 1 data, we make use of the univariate cubic, Hermite interpo- 
lation applied along rays emanating from a vertex and joining to the opposing edge. See 
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Figure 20.44: The data for C l interpolants position and derivative boundary values. 


Figure 20.41. Cubic Hermite interpolation will match position and derivatives at the two 
ends of the interval. We assume that position and derivative information is available on the 
entire boundary of a triangle 7^. 


Si[F](p) = (3 - 2b t )F(V t ) + b?(b t - 1 )F\Vi) 

+(1 - 6i) 2 (26< + l)F(Si) + 6,(1 - b t )-F'(S,) 


where F'(V,) = { ~ — + gllMhl and 

1 — bi 


F'(S, 


. x _ ( x - Xi)Fx(Si) + (y - y^F^S,-) 
*/ ~ I • 


1-fc 


i [F] has the property that it interpolates to 


the boundary data provided by F at V, and on the entire opposing edge e* j . It also matches 
first-order derivatives on this edge and at V,. It does not necessarily interpolate F or its 
derivatives on the other two edges. In order to have an interpolant for the entire boundary 
of the triangular domain, we could try to construct one using the ideas of Boolean sums as 
was done earlier for the C°, side-vertex interpolant. Even though the interpolants Si,Sj , 
and Sk commute so that their Boolean sums are well defined, this approach does not work 
(see [ 1 77]) and so the use of convex combination techniques has been suggested. This leads 
to the interpolant 


or Fl _ W,[F] + b?blSj[F] + bjb]Sk{F} 

1 J - b^ + b]bl+bM 

which has the property that it matches F and its first order derivatives on the entire 
boundary of the triangular domain. In the case where the boundary information has been 
discretized with cubically varying (Hermite) position values and linearly varying cross- 
boundary derivatives, it is possible to obtain a final interpolant with simpler weights in the 
convex combination. Namely, 

slF] = bjb k Si[F] + bjb k Sj[F] + 6jfrS fe [F] 
b%bj + bjbk + b{bk 
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20.3 Tetrahedrizations 

In this section we follow the outline of the previous section as well as possible. Since the 
dimension is one less and since bivariate problems have been considered for a much longer 
period of time, the development in the 3D domain is not as rich as it is in the 2D domain. 
So we cannot parallel the previous section exactly, but most everything generalizes or leads 
to something interesting and often useful. 

20.3.1 Basics 

Definitions, Data Structures, and Formulas for Tetrahedrizations 

Our definition of a tetrahedrization follows very closely to that given for a triangulation at 
the beginning of Section 20.2.1. We start with a collection of points p, — (£,, p», z,), t = 
1, . . . ,N which we assume are not collectively coplanar. We denote this collection of 
points by P. A tetrahedrization consists of a list of 4-tuples which we denote by I t . Each 
4-tuple, i jkl £ I t denotes a single tetrahedron with the four vertices p,- , pj , Pk , and pi The 
following conditions must hold: 

i) No tetrahedron Tijki, ijkl £ It is degenerate. That is, if ijkl £ It then p, , p } , pk , 
and pi are not coplanar. 

ii) The interiors of any two triangles do not intersect. That is, if ijkl £ I t and ap-<5 £ I t 
then Int(7^jA;()r> Int(T a p-,i) = 0. 

iii) The boundary of two triangles can intersect only at a common triangular face. 

iv) The union of all the triangles is the domain D = U^g/, Tijfc/- 

We should point out that condition iii) must hold in the strictest sense, and so tetrahedra 
joining as shown on the right side of Figure 20.45 are not allowed. The reason for this 
condition (and all the others) is the same as the reason for the conditions of a triangulation; 
that is, we eventually wish to be able to define C° functions in a piecewise manner over 
the domain consisting of the union of all tetrahedra. The triangular grid data structure for 
representing triangulations (illustrated in Figure 20.8) generalizes very nicely to a structure 
for representing tetrahedrizations. For example, in Figure 20.46, we show a tetrahedrization 
of the cube into five tetrahedra. We saw earlier in the case of triangulations that once the 
boundary is specified, the number of triangles comprising the triangulation was fixed, and 
moreover, we had a simple approach for determining a formula for the number of triangles 
that existed in the triangulation. This property allowed for the definition of the vectors 
of angles which lead to the criterion for optimal triangulations and was therefore rather 
important. It would be nice if everything extended to 3D in a straightforward manner. That 
is, we would like to say that any polyhedron can be decomposed into tetrahedra and that 
there is a fixed formula of the following form N t = a N b + bN, + c, whereas before, N b and 
Ni are the number of vertices on the boundary and interior, respectively. Unfortunately, this 
is not the case and, in fact, the situation is much worse than that. We saw earlier that any 
polygon-bounded region can be triangulated using only the vertices of the polygon. This is 
one of the first areas where matters differ significantly when going from 2D to 3D. It turns 
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Figure 20.45: The configuration indicated by the diagram on the left is acceptable, while 
that on the right is not acceptable for a tetrahedrization. It is eliminated by condition iii) 
listed above. 
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Figure 20.47: The twisted prism of Schoenhardt [221], which cannot be tetrahedrized. 


out that not every polyhedron can be tetrahedrized. The example illustrated in Figure 20.47 
is originally due to Schoenhardt [221], It can be visualized as a prism which has been 
twisted until each face (a quadrilateral comprised of two triangles) has “buckled” inward. 
Any tetrahedron we form from these vertices must include an edge which lies outside the 
domain of the “twisted prism,” so it is clear that the object cannot be tetrahedrized. 

One very basic operation does carry over in a straightforward manner from 2D to 3D— 
the process of inserting an additional vertex into the interior of an existing tetrahednzation. 
See Figure 20.48. If the new vertex p lies interior to an existing tetrahedron, say T a bcd, then 
this tetrahedron is simply replaced with the four tetrahedron, T a bc P , T a bd P , T b cd P , T ac d P , 
adding a net increase of three tetrahedra. If the new vertex p lies on the common triangular 
face of two tetrahedra, then these two tetrahedra are replaced with six new tetrahedra, 
... rp -]• rp 7 ' 'f resulting in a net increase of four new tetrahedra. 

Tabcp, T bcdp , Ubdp, d aecp , J-ecd P a a edp, resulting m a 

This latter aspect of the number of tetrahedra increasing which is different here from the 
2D case is that net increase in the number of tetrahedra depends on the actual location o 
the interior point to be inserted. This observation points out that not only can the number 
of ways that a data set is tetrahedrized vary, but also the number of tetrahedra can vary. We 
will illustrate this further with some examples that do not even have interior points. 

We have already seen (Figure 20.46) the decomposition of a cube into five tetrahe- 
dra. It is also possible to tetrahedrize the cube into six tetrahedra. This is illustrated in 
Figure 20.49. 

It is interesting to note that from the exterior, the tetrahedrization of Figure 20.49 looks 
exactly the same as that of Figure 20.46 because all external edges are the same. Another 
interesting connection between these two tetrahedrizations of the cube is that one can be 
obtained from the other by “swapping” operations, similar to those used in the Lawson 
algorithm for computing optimal triangulations. Previously, in the case of tnangulations, 
there was the possibility of two triangulations of a convex quadrilateral. The analogous 
situation in 3D is the tetrahedrization of the region formed by five vertices when two tetra- 
hedra meet at a common triangular face. If the line segment joining the two vertices not on 
the common face intersects the interior of the common face, then, analogous to the convex 
quadrilateral case in 2D, there is the possibility of an alternate tetrahednzation. But what is 
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Figure 20.48: Inserting a point interior to an existing tetrahedrization. On the left, the new 
point is interior to a tetrahedron, and on the right it is on a common face of two tetrahedra. 



Figure 20.49: A tetrahedrization of the cube into six tetrahedra. 
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Figure 20.50: Two different tetrahedrizations of five points. 


really different from the 2D case is that the number of tetrahedra changes from two to three! 
This is illustrated in Figure 20.50. This basic operation was applied to the center and upper, 
back-right tetrahedra of Figure 20.46 to arrive at the tetrahedrization of Figure 20.49. 

Another example worth noting in this context is the case where p, = (i, i ,i ),i 
1, . . . , jV. The (Delaunay) tetrahedrization of the convex hull of this set of points con- 
sists of the tetrahedra with vertices Pi , Pi + i , Pj , and Pj+i, of which there are a total of 
((N — 2)(/V — l))/2 tetrahedra. Bern and Eppstein[16] point out that this example pro- 
vides an upper bound on the number of tetrahedra in a tetrahedrization of an JV-vertex 
polyhedron, and that a lower bound is provided by the fact that any tetrahedrization of a 
simple polyhedron has at least N - 3 tetrahedra. 

Some Special Tetrahedrizations 

Following the pattern established in the earlier sections on triangulations, we first discuss 
tetrahedrizations related to Cartesian grids followed by tetrahedrizations associated with 
curvilinear grids. A 3D Cartesian grid (Figure 20.51) involves three monotonically in- 
creasing sequences, x,- , i = 1 , . . . , N x , yj , j = 1, ■ ■ ■ ,N y and zk,k = 1, ... ,N Z . The 
grid points have coordinates (x, , y 3 , z*) and these points mark out a cellular decomposition 
of the domain consisting of regular parallelepipeds. Each of these cells can be tetrahedrized 
in a manner similar to that given for the cube in the previous section. Probably the most 
popular is the tetrahedrization involving five tetrahedra shown in Figure 20.46. So as to 
not end up with a nontetrahedrization with problems similar to those shown on the right 
side of Figure 20.45, it is necessary to “alternate” the tetrahedrization from one cell to the 
next so that adjoining cells have the same diagonal on the common faces. This alternate 
tetrahedrization is not really different but is just a rotation of its companion. It is shown in 
Figure 20.52. Another popular choice is the tetrahedrization shown in the upper-left comer 
of Figure 20.53. It has the advantage that all of the tetrahedra are the same shape (up to 
mirror images). Actually, it turns out that there are six different tetrahedrizations of a cube 
(parallelepiped). See Nielson [183]. We have previously shown pictures of two of them m 
Figure 20.46 and Figure 20.49. The other four are shown in Figure 20.53. 

All six tetrahedrizations of the cube are comprised of five primitive tetrahedra, which 
are shown in Figure 20.54. We use the names OF, IF, 2Fr, 2F1 and 3F for these tetrahedra to 
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( x j, Yj , z k ) 

Figure 20.51: Three-dimensional Cartesian grid. 
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Figure 20.54: The five primitive tetrahedra comprising the tetrahedrizations of the cube. 


indicate the number of exterior faces for each tetrahedra. There are two different primitive 
tetrahedra with two exterior faces; one is a mirror-image version of the other and so it 
cannot be rotated to the other. The tetrahedron OF has volume 1/3 and all the others, liave 
volume 1 / 6 . During informal discussion we most often use the names 3F = “comer,” 2Fr 

or 2F1 — “right wedge” or “left wedge,” IF — “kite” and OF = equi or fatboy. 

In a joining similar to that shown in Figure 20.50, three 1 F tetrahedra can come together 
to form the same exact shape formed by a OF and a 3F together. Also a 2F1 and 2Fr together 
form the same shape as a IF and a 3F, but two 2Fr’s or two 2Fl’s cannot share a common 
face and remain inside a unit cube. There are four tetrahedrizations (each composed of 
three primitive tetrahedra) of the prism making up half of the cube. They are 3F, IF, 2F1, 
3F, IF, 2Fr; 2Fr, 2F1, 2Fr; and 2F1, 2Fr, 2F1. In Figure 20.55 we show the dual graphs of the 
six tetrahedrizations of the cube. A node is a primitive tetrahedron and an arc is a common 
triangular face. As expected, in each case the names add to twelve. 

Each of these six tetrahedrizations has unique and interesting properties. The tetra- 
hedrization of Figure 20.46 and Figure 20.49 both “swap” diagonals on all three pairs of 
opposing faces. The tetrahedrization shown in the lower right of Figure 20.53 swaps the 
diagonals of two pair of opposing faces and that of the upper right swaps one pair. The 
two tetrahedrizations on the left of Figure 20.53 do not swap any diagonals of any op- 
posing faces. The tetrahedrization of the upper left of Figure 20.53 can be realized with 
three cuts of the entire cube, while the others cannot. This particular tetrahedrization also 
has the unique property of being composed of only 2F primitives whose faces are all right 
triangles, and all six of them share the diagonal of the cube as a common edge. This tetra- 
hedrization has been discussed and used widely. It is called the CFK triangulation of the 
cube after Coxeter [47], Freudenthal [79], and Kuhn [137], A replacement rule can be 
used to generate this tetrahedrization. Using the labeling scheme of Figure 20.46, we start 
with the four vertices P 2 ._ lt i = 0, 1,2, 3 and replace each vertex V jt other than V'o and 
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Figure 20.55: The six tetrahedrizations of the cube shown as dual graphs. (These are the 
only tetrahedrizations of the cube.) 
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Figure 20.56: Face triangulations that are not consistent with any tetrahedrization of the 
cube. 


V 7 , with V j+ i + Vj-i-Vj. Explicitly, this will successively generate the six tetrahedra: 
P0P1P3P7; P0P2P3P7; P0P2P6P7; P0P4P6P7; P0P4P5P7; P0P1P5P7. The CFK tnangula- 
tion generalizes to n-dimensions as does the “replacement” algorithm for generating the 

simplicial decomposition. . 

It is interesting to note that not all possible face triangulations are realized by the six 
possible tetrahedrizations of the cube. In addition to the five different face triangulations 
which are realizable (note that two tetrahedrizations have the same face triangulations) 
there are three others which cannot be realized. They are shown in Figure 20.56. In or- 
der to determine these eight unique face triangulations, we start with the 64 = 2 face 
triangulations and then group them into these eight equivalence classes by rotations. 

Theorem: It is impossible to tetrahednze a cube and yield face triangulations as 

shown in Figure 20.56. 

Proof: We give only the proof for the case in the top center, as the others are similar. 
We use the same labeling as shown in Figure 20.49. We start with the face 457. Only vertex 
0 can be attached to the face 457, which gives the tetrahedron 0457. The internal face 047 
must be shared by some other tetrahedron. Any vertex, however, cannot be joined to the 
face of 457 without violating the conditions of the face triangulations, so this completes the 
argument. 

Earlier we discussed triangulations related to curvilinear grids. We now take up the 
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Figure 20.57: Single cell of a 3D curvilinear grid, 
topic of tetrahedrization ot 3D curvilinear grids. Analogous to the 2D situation, a 3D 

curvilinear grid is specified by three geometry arrays i- l,... j _ 

l,.. . ,Ny- k = 1, . . . , N s . In the 2D case a cell C,j consisted of the quadrilat- 
eral with vertices (x tj ,y tj ), (x i+1J ,y i+lJ ), (x iJ+u y tJ+1 ), (x i+lJ+uyi+1J+1 ), and the 
cells serve as a decomposition of the domain. 

In the 3D case, matters are not as straightforward as we might expect, and there are 
some areas where we need to be concerned. These have mainly to do with just exactly 
what comprises a cell. In 3D the cell C ijk has the eight vertices (x afcc , y abc , z abc ), a = 
i, i + 1, b — j, j + 1, c = k, k + 1, but there is not always a consistent definition for the 
cell boundaries. We mention briefly some possible choices. If the geometry arrays are 
constrained so that each collection of four vertices of the six “faces” of the cells are copla- 
nar, then an obvious choice for the cell boundaries is this common planar quadrilateral. In 
this case the cells are hexahedron, and it is relatively easy to determine whether or not an 
arbitrary point (x, y, z) is in a particular cell or not. Often this planarity condition does 
not hold, and cell boundaries are taken to be the parametrically defined (hyperbolic) sur- 
face obtained by substituting 0 or 1 for any of the parameter values s, t . u in the following 
trilinear mapping: 

Oij.k(s,t,u) = (1 — s)(l - t)(l-u)P iJ:k + (1 -s)(l -t)uPi ’ iJtk+1 

+(1 - s)t(l - u)P it j + i k + (I - s)tuPi j +l k+l 
+s(l - t)(l - u)P i+lJ k + S (1 - <)«P, +Ijifc+1 
+St(l - u)P i+ l,j + l,fc + stuP i+li j + l k+l 

where 

P >J,k = ( x i,j,k,yiJ,k,Zi,j,k) 

Given a point (x,y,z) in the cell Cijk, the value (s, t, u) which associates with it via 
the trilinear mapping is called the corresponding computational coordinate. In fact, in 
order to determine whether or not an arbitrary point is in this type of cell or not requires 
that we solve the three nonlinear equations which represent this association. This can 
be a considerable problem from a computational point of view. Most methods use some 
heuristics to obtain an initial approximation for some type of Newton’s method. Another 
choice for the cell boundaries, in the event the four vertices of a face are not coplanar, 
is to choose them to be piecewise planar. That is, a diagonal edge is selected and the 
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Figure 20.58: A curvilinear grid cell (polyhedron) that cannot be tetrahedrized. 


boundary between the two cells consists of the two triangles which result. Often the cell 
would be further decomposed into tetrahedra, thus leading to an overall tetrahedrization of 
the curvilinear grid. We should point out that not all choices for the diagonals can lead to 
a tetrahedrization of the cell. In order to be specific about this, consider the cell illustrated 
in Figure 20.58. This cell was created from a unit cube by cutting notches in the faces 
so as to force the diagonal edges P2P7,P4Pi>P3P5,P3POiPoP6,P6Ps to be exterior to the 
cell If the depth of the notches is e, then this results in the points p 0 = (0, e, 0) , p i = 
(1 — e, 0, e), P 2 = (e, 1, e), p 3 = (1, 1 - e, 0), p 4 = (e, 0, 1 - e), Ps = (M. 1), Pe = 
(0, 1 - e, 1), P 7 = (1 - e, 1, 1 - e). Note that p 6 ,P3,P4 and pi all lie in the plane 

x -j- z — 1=0 and P 2 , P 7 , Po and p$ are in the plane x — z — 0. 


Theorem: The polyhedron of Figure 20.58 cannot be tetrahedrized. 

Proof: Consider the triangle face with vertices pe , P4, and pi- In any tetrahedrization, 

this face must be joined to some vertex to form a tetrahedron. By considering the remaining 
five vertices p 5 , p 0 , p 2 , Pi , and p 3 we find that only p 3 would not lead to a tetrahedron with 
an edge which is outside the cell. If the tetrahedron p 6 , P4, P7, and p 3 is included in the list 
of tetrahedra, then the interior triangle face p 3 p4P7 must connect to another vertex (besides 
p 6 ) to form a tetrahedron. But a consideration of each of the possible vertices Ps ■ Pi , P 2 , 
and po each lead to an edge which is exterior to the cell and this concludes the argument. 

We conclude this discussion on the tetrahedrization of the cells of a curvilinear grid 
by pointing out that some hexahedra will decompose into seven tetrahedra. Consider the 
cell of Figure 20.57 and let the six faces be planar, but assume that the four diagonal 
points p ijk , p.+u+m+i, PijM i andp t+u+ i,* are not coplanar, so that they will form 
a tetrahedron. Remove this tetrahedron, leaving two prisms with two planar quadrilateral 
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faces which can each be decomposed into three tetrahedra. We should point out that we 
have observed cases where this decomposition was the Delaunay tetrahedrization. 

In Section 20.2.1 we described two different approaches leading to nested subdivi- 
sion triangulations and pointed out their potential value in multiresolution approximations. 
These both have analogs in 3D and are shown in Figures 20.59 and 20.60, respectively. 
The first one is based upon recursive subdivision and the second one is called “symmetric” 
subdivision and is related to the CFK tetrahedrization of the cube [170]. It is composed of 
six 2Fr’s and two 2FFs and is the same shape and twice the size of one 2Fr. 

It should be noted that if primitive tetrahedra of the shape shown in Figure 20.6 1 are as- 
sembled as in Figure 20.60, then we obtain a composite tetrahedron which is twice the size 
and exactly the same shape as the primitive tetrahedron. This particular tetrahedrization of 
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(0, a, a/2) 



(a, a, 0) 


(a, 0, a/2) 


Figure 20.61: A tetrahedron that can be tetrahedrized into eight tetrahedra, each of which 
are the same shape as the original yet half the size. 



needle ca P 



wedge sliver 


Figure 20.62: Examples of poorly shaped tetrahedra. 


tetrahedra is related to the Delaunay tetrahedrization of the BCC lattice, which is the union 
of the lattices {(i, j, k) : i, j, and k are integers) and {(i + j + \ , k + i) : *, j, and k are 
integers). See also Senechal [229] for a discussion of tetrahedra that can be decomposed 
into similar tetrahedra. 


20.3.2 Algorithms for Delaunay Tetrahedrizations 

Analogous to the examples of Figure 20.19, examples of poorly shaped tetrahedra are 
shown in Figure 20.62. The sliver has small dihedral angles, but need not have any small 
planar angles. Several measures of the quality of tetrahedrizations have been proposed. See 
Baler [12] and Field [86]. Take as an example the ratio of the inradius (radius of inscribed 
sphere) and the circumradius. The problem here is that there is no apparent way to order 
the collection of all tetrahedrizations of a point set. The approach of lexicographically or- 
dering the associated vectors of angles, as we described in Section 20.2.2, does not extend 
to 3D because the number of tetrahedra in a tetrahedrization is not necessarily fixed. Nev- 
ertheless, the Delaunay tetrahedrization of the convex hull which is dual to the Dirichlet 


480 


Tools for Triangulations and Tetrahedrizations 



Figure 20.63: Different cases of swapping for 3D version of Lawson’s algorithm. 


tessellation is well defined (in the absence of neutral cases where points lie on a common 
sphere), so the remainder of this section is devoted to a discussion of the extension of the 
previously discussed 2D algorithms for computing the Delaunay tri angulations to the case 
of 3D tetrahedrizations. 

Extension of Lawson’s Algorithm (Incremental Flipping): It is possible to ex- 
tend this algorithm to 3D, but the extension is not as simple as one might expect. The first 
major difference that one encounters is the character of the basic swapping step. In 2D we 
take an edge and consider the quadrilateral formed by the two triangles which share this 
edge. If the quadrilateral is convex, we can swap the diagonal if this step moves us closer to 
the optimal solution, which can easily be determined by applying the circle inclusion test. 
Two triangles are replaced by two other triangles. But the analogous steps in 3D can lead 
to a situation where the two tetrahedra sharing a face can be replaced w ith three tetrahedra. 
See Figure 20.63 for an example. 

Joe [122] showed that if the points are inserted in a particular manner, then incremental 
flipping will lead to the optimal Delaunay tetrahedrization. Edelsbrunner and Shah have 
generalized these results [72]. See also [82]. Software based upon these ideas is provided 
by the Software Development Group at the National Center for Supercomputing Applica- 
tions and is available at the World Wide Web site: 

http://www.ncsa.uiuc.edu/SDG/Brochure/Overview/ALVIS.overview.html. 

Extensions of the Algorithm of Green and Sibson: There does not seem to be 
an apparent method of extending this type of algorithm to 3D. The algorithm is dependent 
upon the “contiguity list,” and here lies the difficulty to extend to 3D. We included this 
algorithm in our selection of 2D algorithms so that this very point could be made. Some 
concepts extend easily to 3D and others do not. 

Bowyer’s Algorithm for 3D: It is a straightforward exercise to extend Bowyer’s 2D 
algorithm to 3D. In fact, the original paper of Bowyer [21] describes the algorithm for 
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arbitrary dimensions. Bowyer also mentions that with some care, the algorithm can be 
extended to other domains. In [164] there is a brief discussion of Bowyer’s algorithm 
along with some code. 


Watson’s Algorithm for 3D: The original description of Watson’s algorithm applies 
to arbitrary dimension. In Watson’s paper [254] results for 2, 3, and 4 dimensions are 
reported. Information on implementing this algorithm in 3D is given by Field in [86] and 
[87], It is also the basis for the 3D algorithms discussed in [29]. 

Embedding/Lifting Algorithms for 3D: Software for computing general dimension 

convex hulls and Delaunay tetrahedrizations, based on the relationship mentioned earlier 
in Section 20.2.2, are provided by the Geometry Center, University of Minnesota at the 
WWW site: 

http://freeabel.geom.umn.edu/software/download/qhull.html. 

20.3.3 Visibility Sorting of Tetrahedra 

We first give a motivation for the definition of and the need for a visibility sort. We use 
the example of volume rendering, which is a means of graphing (visualizing) a density 
function (cloud) S(x, y, z) defined over a 3D domain (which is often a cube). A viewpoint 
V is selected along with a projection plane. A rectangular portion of the projection plane 
is subdivided into a rectangular array of subrectangles which associate directly with the 
pixels of an image to be generated. The value for each pixel is defined by 

F(i,j) = f° 6(s)C(s)e~ ^ ds + F 0 e- ‘M* (20.3) 

J o 

where the integral is taken along the ray emanating from the viewpoint and passing through 
the center of the subrectangle associated with the pixel at location ( i,j),F 0 is the back- 
ground intensity, and D is a distance along the ray sufficiently large so that the ray com- 
pletely passes through the domain of interest. The function C, also defined over the same 
domain as J, is called the color function and governs the color of light emanating (by re- 
flection, say) from a point within the density cloud. In actual application the integrals are 
approximated by numerical schemes based upon sampled values of the integrand. The 
sample values are often obtained by some simple interpolation into the cells covering the 
domain And these cells are often a result of the positions where S has been measured. If 
we let 0 = xo < *i < * 2 < ■ • ■ < *»-i < x n = D be the distances from the viewpoint 
to each sampled value along the ray, then the upper Riemann sum approximation to this 
integral is 

n (20.4) 

i=0 j=« + 1 

where C{ = C(xi) y tj = e ~ Ax * 5 ( x ^ and Ax* = x t - - x*_i . This discrete approximation 
can be computed by the compositing process 

Fi — riTi-i + hy 


(20.5) 
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Figure 20.64: An example of three tetrahedra that cannot be visibility ordered. 


where U — A 

Another way to view this compositing process is as a simple model of transparency, 
where an object of thickness Ax t attenuates the incoming light intensity F, _ t by the factor 
U, and this object emits light of intensity /,. Algorithms which accumulate these values 
into a frame buffer (with each location holding the value for a pixel) can either be image- 
space-oriented or object-space-oriented. Image-space algorithms proceed along the lines 
of our development here and accumulate all contributions for a pixel along a particular ray. 
Object-space algorithms compute exactly the same values but the calculations are done in 
a different order. These algorithms sequentially process each cell by accumulating into the 
proper location of the frame buffer all contributions of a particular cell. Due to the nature 
of the compositing process, it is mandatory that these accumulations be done in the proper 
order. It is this latter approach which motivates the definition of and need for visibility 
sorting in this context. 

Definition of Visibility Order: Let T and V be tetrahedra of a tetrahedrization and 
let V be the center of perspective projection. If there is a ray emanating from V which 
intersects T ( before T, then T is said to precede T f and we write T <T* . 

The purpose of a visibility sort is to find a linear ordering of all of the tetrahedra of a 
tetrahedrization so that the ordering relation is never violated. 

Definition of Visibility Ordering: A visibility ordering of a tetrahedrization is a se- 
quence, rix , n 2 , . . . , 7 it which has the property that whenever T n% < T Uj then i < j. The 
implication of the definition of visibility ordering for splatting or object-space traversal 
algorithms for volume rendering is that a tetrahedron T must be processed (sampled and 
composited into the frame buffer) before T* whenever T < T f . 

A couple of items should be noted at this point. The relation of visibility order is not, 
in the strict mathematical sense, a partial ordering. A partial ordering is required to be i) 
transitive: x < y, y < z implies x < z; ii) antisymmetric: x < y and y < x implies 
x = y; and iii) reflexive: x < x. It is entirely possible that a visibility order could not exist 
at all due to the presence of cycles as shown in Figure 20.64. 

Knuth [136] has discussed in some detail (including MIX programs) the topological sort 
algorithm as a means of “embedding a partial order in a linear order.” A linear ordering 
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Figure 20.65: An example of the topological sort algorithm. 


is a partial ordering where either x < y or y < x for all x , y. Even though this does not 
strictly apply in the context of a general tetrahedrization, the basic ideas (mainly due to the 
manner in which it is described) are very useful for developing visibility-sorting algorithms 
for specific applications, so we include a description of the topological sort algorithm here. 

Topological Sort Algorithm: The topological sort algorithm as described by Knuth 

[136] starts with a directed, acyclic graph (DAG). The DAG can be represented with a 
diagram using nodes and arrows. See Figure 20.65. The nodes represent the elements of 
the set to be ordered, and an arrow from node x to node y represents the relation of the 
partial ordering, x < y. The algorithm is simple. Any node that has no incoming arrow is 
removed from the DAG (with all of its attached arrows) and placed in the linear ordering. 
This process is repeated until the DAG is empty. It is easy to prove (left to the reader) that 
if the DAG represents a partial ordering, a linear ordering will always be produced by this 

algorithm. ... , 

Max [166] has discussed the application of the ideas of the topological sort algorithm 

to the problem of producing a visibility sort for a cellular decomposition of a domain. Max 
defines the order relation in the following way. The DAG contains an arrow for each face 
common to two cells x and y. The arrow is directed from x to y if the viewpoint is on the 
same side of the face as x, meaning that y must be processed before x. Max mentions that 
the topological sorting algorithm will be successful “if every ray through the data volume 
intersects it in a single sequence of adjacent cells.” Of course, if the cell complex contains 
cycles (see Figure 20.64), then a visibility sort is not possible. Williams [257] discusses 
similar algorithms applied to a very general cellular decomposition which may contain 
empty cavities. 

We conclude this section with some rather interesting properties about die special case 
of the Delaunay tetrahedrization of the convex hull of a collection of 3D points. The power 
of a tetrahedra is defined as D 2 - R\ where D is the distance to the viewpoint from the 



484 


Tools for Triangulations and Tetrahedrizations 



Figure 20.66: Elements of the definition of the power of a tetrahedron. 


center of the circumsphere of the tetrahedron and R is the radius of the circumscribing 
sphere. See Figure 20.66. A visibility sort can be accomplished by a simple sort based 
upon the power. This property is covered in [69] and used by Max, Hanrahan, and Crawfis 
[167]. We caution the reader that this approach breaks down in the presence of neutral cases 
where possibly several tetrahedra have the same power (as in the case of decomposing the 
cube). One additional interesting observation in this context is that a sort based upon the 
power of the tetrahedra does not require the neighborhood information that is required for 
the algorithms using the ideas of topological sorting. Another method which does not use 
adjacency information is described by Stein, Becker, and Max [240], 

20.3.4 Data-Dependent Tetrahedrizations 

Lee [148] has investigated the topic of data-dependent tetrahedrizations. This work gener- 
alizes from 2D to 3D the ideas and techniques of [67] and [225], Similar to the algorithms 
of [225], simulated annealing is used. The initial tetrahedrization is the Delaunay tetra- 
hedrization of the convex hull of the independent data site locations. Local swapping of 
tetrahedra is performed based upon random values compared to an annealing schedule and 
a cost function. This ‘‘randomness” of the simulated annealing approach allows the algo- 
rithm to escape local extrema of the cost function. Local swapping for 2D simply involves 
the choice of one or the other of the diagonals of a quadrilateral. In 3D the situation is more 
complex. There are four cases shown in Figure 20.63 which are the same as those used in 
the 3D version of Lawson’s algorithm. In the first case, three triangles are swapped for two. 
The second case is the reverse of the first — two tetrahedra are replaced by three. The third 
case is where two triangles are on the boundary of the convex hull and the two tetrahedra 
can be swapped for two other tetrahedra. In the last case four tetrahedra are swapped for 
four other tetrahedra. 

In Section 20.2.4, we described the cost function used by Dyn, Levin, and Rippa [67]. 
Analogous to these cost functions for 2D, Lee [148] uses the following criterion for 3D: 

Gradient Difference: Let T\ and T 2 be two tetrahedra with a common triangular face. 
Let G i be the gradient of the linear function which interpolates the data at the four vertices 
of T\, and let Go be the similar gradient for the linear interpolant of T 2 . The gradient 
difference is defined as ||Gi — G r 2 ||- 
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Method 

RMS Error 

Delaunay 

.007475 

Difference in Gradient 

.005445 

Jump in Normal Derivative 

.004361 


Figure 20.67: Errors for the piecewise linear interpolant using different tetrahedrization. 


Jump in Normal Direction Derivatives: Let L\{x, y, z) — aix4&iy4ci*4di be the 
linear function which interpolates to the data at the four vertices of T x andletL 2 (z, y,*) = 
clox 4 b 2 y 4 c 2 z 4 d 2 be the similar function for T 2 • Let N — {ti x , n y , be the normal 

(normalized) of the common triangular face of 7\ and T 2 . T>i = ain* 4 &in y 4 c x rt z 
is the directional derivative of L x in the direction of N. D 2 = a 2 n x 4- b 2 n y 4 c. 2 n z 
is the analogous value for To. The jump in normal direction criterion is \D X — D 2 1 — 
|(ai — a 2 )n x 4 (61 — b 2 )n y 4 (ci — c 2 )n z )\. 

Some example results reported by Lee [148] are repeated here in Figure 20.67. This 
example involves a test function, F(z, y t z) = (tanh(9y - 9x - 9 z) 4 l)/9, which provides 
the dependent data. The piecewise linear interpolant over the tetrahedrization is compared 
to the test function. The RMS errors are based upon evaluations of the functions and this 
approximation over a 20 x 20 x 20 Cartesian grid. The dependent data site locations are 
taken to be 1000 random points in the unit cube. 

In Figure 20.68 are some graphs which can be considered as 3D analogs of the graphs 
shown in Figure 20.31 of Section 20.2.4. Similar to the 2D case, the data-dependent tetra- 
hedrization involves some badly shaped tetrahedra. This is the cost of having an optimal 
(or nearly optimal) piecewise linear approximation. 


20.3.5 Affine Invariant Tetrahedrizations 


In this section we extend the results of Section 20.2.5 on affine invariant triangulations 
to that of affine invariant tetrahedrizations. Prior to discussing the characterization and 
computation of this type of tetrahedrization, we make some comments about the need for 
such a tetrahedrization over and above those reasons for the 2D case. It appears that as the 
dimension of the independent data increases, our need to be concerned about lack of affine 


invariance also increases. 

One source of 3D independent data is the case of time-varying 2D data. In some cases 
the data measurement locations might stay fixed over time and in some cases they may 
vary over time. Let us say, for example, that we have a vector field which is known (by 
means of a numerical simulation) at the locations of a 2D curvilinear grid ( Xij , Vij), * — 
l t N x ; j = 1, . . . , N y . As time proceeds, the vector field varies, but the dependent 
data site locations stay fixed. So in this case, we have data which can be represented as 
{Vijk',Xij,yij,t k ), i=l,...,N x , j = 1,... ,N y , k = If the definition of 

a modeling function F(x } y,t), designed to interpolate the data, Fyxij.yijytk) — fijk, 
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Figure 20.68: Data-dependent tetrahedrization compared to the Delaunay tetrahedrization. 
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is based upon a tetrahedrization of the 3D independent data {xij , yij , t/c), then this model 
will not necessarily be affine invariant, and the units used to measure and represent the 
physical coordinates and time could have an effect on the modeling function F{x , y, t) 
and subsequently an effect on the visualization and analysis. The same problem could 
also occur for a time-varying vector field over a curvilinear grid which also varies over 
time — that is, data of the type (Fy/t; Zijfc, yyfe.ffe). » = 1, • • • . Nr, j = 1, ■ • ■ , F y , k = 
1, . . . , JVf. In general, any tetrahedrization of the independent data of , yj , ) -» 

where the choice of the units of measurement used for the independent data could lead 
to a nonuniform scaling, could have the problem of being dependent on the choice of the 
units used. If each of the variables use the same units there will be no problems of this 
type, because a scale transformation of the form x <- ax, y <- ay, z f- a z, where the 
scale change is uniform for each variable, will not affect the tetrahedrization. It is only the 
nonuniform scaling x ax, y ^ by, z ^ cz which creates the problem. An example 
of a scale change affecting the tetrahedrization is shown in Figure 20.69. Here there are 1 0 
data points. In the right image, the data has been scaled in the y variable by a factor of 2. 
Not only does the tetrahedrization change, but even the number of tetrahedra changes. The 
Delaunay tetrahedrization of the original 10 data points has 18 tetrahedra and the scaled 
data has 13 tetrahedra. 

We now describe the 3D version of the affine invariant norm, which leads (by way of 
the Dirichlet tessellation) to an affine invariant tetrahedrization. Actually, we can define it 
so that it is clear what the generalization is for any dimension. Let 

\\{x,y,z)\\l = {x, y ,z)(vvr l ^ y 

where V is the 3 x N matrix of translated data values 

( Xi~Hx X2 — fix ■■■ XN — fix \ 

t/i -fiy Vi-Vy • • • Wf - A*» I • 

Zl - fi z z 2 - flz ■■■ Z N - fiz J 

As with the 2D case, there are some different approaches to modifying an existing tetra- 
hedrization procedure. Probably the simplest is to preprocess the data with the transforma- 
tion given by the lower triangular matrix, L(V) which results from the Cholesky decom- 
position of (VV*) -1 

L(V)L(Vy = (VV')- 1 . 

Explicitly in the 3D case, we use the transformed data 

Xi = h\Xi +hiVi + hi*> 

Yj = hlVi + ^32 z i 
Zi = 


— _ a n a 13 

U1 — V a iu *21—1 ’ *31 / 1 

/ll *11 


where 
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Figure 20.69: Delaunay tetrahedrization of 10 data points and a scaled version of the same 
data points. 
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122 = \/«22 - (^21 ) 2 > ^32 = 


Q 33 ~ ^21^31 
/ 22 


/33 = \/ a 33 “ (^3l) 2 — (^32 ) * 

( aij ) = (W)~ l = L(V)L(V)* 


P 2 P 2 _ V 2 

1 / yl x 

— — I — (E ry E~ — E r jEyz) 
aet \ v v V- 

Z^xy^yz ^xz^y 


-(E 


ary 


E J 


.SvO 


V V 
^ry ^yz 


V E 2 

^ I2 y 


V2V2 __ ^2 


fV V2 _ V V 

\^yz2-J X <-*xz '- t xy ) 


(Ey Z ^ x — E r2 E ry ) 




■ 5> 
■ry 


where 


£* 

S y 

E? 


N 

N 

eL(*.-a0 2 

jv 

TT'W ( _ .. W, 




Sry “ 


Enz — 


= 


iV 


.iv 


AT 


^ N 


Vx = 
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V ' tV 

Ei=i w 

A*y = 

N 


sr N 
L*i= 1 ** 

P* = 

N 

-fly) 
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-A**) 


-Pz) 



N 


and 


det = E 2 E 2 E 2 + 2Ej 


^lyEyzEri - E^(Eyj) 2 - S 2 (Eu) 2 E 2 (E a :y) . 

We conclude this section with some examples illustrating this affine invariant norm and 
its use in characterizing affine invariant tetrahedrizations. In Figure 20.70 there are four 
graphs of 1 3 data points. The transparent ellipsoids represent all the points that are 0.50 and 
1 .0 units from the center point using the affine invariant norm. The different graphs show 
the data after it has undergone an affine transformation. The original data is displayed in the 
upper left. The upper right shows the data after it has been rotated by 44 degrees about the 
2 axis. The lower right is after it has subsequently been scaled in the x variable by a factor 
of 1.5. The lower left is after it has been scaled in y by a factor of 0.6. A close examination 
of these graphs will show that the relative distances (as measured by the affine invariant 
norm) between points is unchanged by these transformations. Figure 20.71 shows an affine 
invariant tetrahedrization. In comparison, the conventional Delaunay tetrahednzation is 
shown in Figure 20.72. 


20.3.6 Interpolation in Tetrahedra 

As with the bivariate case covered in Section 20.2.6, there are two concepts of interest for 
interpolation in tetrahedra. The first is concerned with the amount of boundary data that 
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Figure 20.70: Examples illustrating the affine invariant norm. The ellipsoids are 0.50 and 
1 .0 units from the center point. 




- n . no 


Figure 20.71: Examples of affine invariant tetrahedrization. 





Figure 20 . 72 : Delaunay tetrahedrization of the same data as in Figure 20.71 . 
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is provided or assumed to be available. This can be discrete data provided at a finite num- 
ber of locations (usually the vertices or midpoints) or transfinite data where boundary data 
values are assumed to be available at all locations on the boundary. The second concept 
relates to the degree of continuity of a piecewise defined interpolant using the local inter- 
polants described here. C° interpolants use only boundary position data and lead to overall 
interpolants which are continuous. C 1 interpolants utilize first-order derivative information 
and lead to global interpolants which have all first-order derivative continuous. These two 
concepts lead to four possibilities which we discuss below. 

C°, Discrete Interpolation in Tetrahedra 

Analogous to the bivariate linear interpolant which will match predescribed values at the 
three vertices of a triangle, there is a unique trivariate linear interpolant which will match 
data at the four vertices of a tetrahedra, T{j k i. Given F (Vi) , F (Vj) , F (V) c) and F(V5), the 
coefficients of this linear function which interpolates this data 

F(x } y, z) = a + bx + cy + dz 
can be found by solving the linear system of equations 


a 4 bxi 4 cyi 4 dzi 

= F{Vi) 

a 4 bxj 4 cyj 4 dzj 

= F(Vj) 

a 4 bxk 4 cyk 4 dz k 

= F{V k ) 

a 4 bxi 4 cyi 4 dzi 

= F(V,) 


As before, it is also possible to use barycentric coordinates. The barycentric coordinates of 
a point V = (x, y, z) are defined by the relationships 

V — bi Vi + bj Vj 4- bk Vk 4- bi V[ 

1 = bi + bj bk + bi 

and the linear interpolant has the form 

F(x, y, z) = F(V) = biF(Vi) 4- bjF{Vj) 4 b k F(V k ) 4 b ,F(Vi). (20.6) 

As before, there are several ways of defining or computing barycentric coordinates. The 
analog of the ratios of areas we saw before is the ratio of volumes of subtetrahedra, 

L V ol(Tv jki) u _ Vol(TiVkt) i _ V ol{Tjjvi) , __ Vol(Tjj k v) 

{ ~ VolWjkt)' bj ~ Vol(Tijki) ’ * " Vol(T ijk i) ’ j ~ Vol(T ijkl ) 

where Tyjki is the tetrahedron with vertices V, Vj } 14, and V} and similar definitions for 
the other subtetrahedra. The volume of a tetrahedron, T a & C d, with vertices a, b t c, and d is 

Vol(T abcd ) = i [{d - o) • ((6 - a) x (c - a))] 




494 


Tools for Triangulations and Tetrahedrizations 


V, 



v i 


Figure 20.73: Data site locations for trivariate quadratic interpolation. 


Also determinants can be used, 
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Given the values at the four vertices and the six midpoints of a tetrahedron, there is a 
unique trivariate quadratic which interpolates this data. 


Q{x, y,z) = F(Vi)bi(bi - bj - b k - 6/) + F{Vj)bj{bj - b x - b k - b t ) 
+F(V k )b k (b k - b z - bj - b t ) + FiV^btibt - b t - bj - b k ) 
+F{M lk )4b t b k + F(Mji)4bjbi + F(A/y)46 i 6 i (20.7) 

+F(M jk )4bjb k + F(Mu)4bib l + F{\f ki )4b k b t 

where M zj — {Vi 4- V})/2 and the other midpoints are defined similarly. See Figure 20.73. 


C°, Transfinite Interpolation in Tetrahedra 

As before in Section 20.2.6, we give a sampling of interpolants. One is a generalization of 
the side- vertex interpolant and the other is a generalization of the C* interpolant. Both of 
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Figure 20.74: Notation for the face-vertex interpolant. 


these bivariate interpolants were discussed previously in Section 20.2.6. 


The C°, Face- Vertex Interpolant: Analogous to the basic interpolants used to con- 
struct the side-vertex interpolant, we have the interpolants which consist of linear interpo- 
lation along edges joining a vertex and the opposing face 


At[F\ = biF(Vi) + (1 - 6,)W) 

Aj[F] = b j F(Vj) + (l-b j )F(F j ) 

A k [F] = b k F(V k ) + (\-b k )F(F k ) (20.8) 

A,[F] = biF(Vi) + (1 — bi)F(Fi) 


where F{ = 


bj Vj+bkVk+bjVi 

bj + bk 4- bi 


Fj = 


bi 


, „ biVi+bjVj + b.Vk 

andF,aS 6, + bi , + b k — 
four interpolants leads to 


See Figure 


Vi+b k V k +b,V, _ 6, Vi + bjVj +biVi 

bi + b k + bi ' k bi + bj + bi 

20.74. Computing the Boolean sum of these 


A[F] = (1 - bi)F(Fi) + (1 — bj)F(Fj) 4- (1 — b k )F(F k ) + (1 —bi)F(Fi) 

— (b k + bi)F(S k i) — (bi + bi)F(Su) — (bj + bi)F(Sjt) 

-(bj + b k )F(S jk ) - (bi + b k )F(S ik ) - (bi + bj)F(Sij) (20.9) 

+biF(Vi) + bjF(Vj) + b k F(V k ) + b,F[Vi) 

where + b ' - - , mn = kl,il,jljk,ik,ij. 

b m + v n 
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The C* Interpolant (for a tetrahedron): The analog of the bivariate C* interpolant 

described in Section 20.2.6 is 


C'[F] = biFiVD+bjFW + bkFW + btFiVt) (20.10) 

+W | F(Q t ) - F(Vi) - (b, + F(V : ) - [b k + F(V k ) | 

+W k |f(Q*) - (ft,- + F(K) - (ft; + y) F(VS-) - (fti + y) F(V})j 

+«-; { F(Qj) ~ (bi + F(Vi) - (ft, + F(I4) - (ft, + F(V)) j 

+Wi {f^) - (ft; + F(l$) - (ft, + F(V h ) - (ft, + F(Vt ) | 

where Q, = (ft, + y) V, + (ft; + j-) Vj + ( b k + — ) V k , 


W, = 


2 / 6, bj b k 

(36, -f 6/)(36j? + bi)(3bft + 6/) 


manner. 


and the other Q's and W's are defined in a similar 


C 1 , Transfinite Interpolation in Tetrahedra 


The C 1 , Face- Vertex Interpolant: It is a straightforward process to extend the C 1 , 

transfinite side-vertex interpolant to a tetrahedral domain, Tijki. It is called the C 1 , face- 
vertex interpolant and we assume that position and derivative information is available at all 
locations on the four faces which make up the boundary of the tetrahedron T^ki . The basic 
face- vertex operator is defined as 

S,[F](p) = ft, 2 ( 3 - 2 bi)F(Vi) + ft? (ft, - l)F'(Vi) 

+(1 -i,) 2 (26 t + l)F(Fi) + &i(l - ft,-) 2 F'(F) (20.11) 


where F'K) = (^ - ^(^) + fa - + ( ^ ^ 

1 - bi 


F‘(Fi) = 


(g - xj)F x (Sj) + (y- yi)Fy(Sj) + (z - z i )F z (S l ) 

1 — bi 


The point F, is the inter- 


section point of the ray from Vi through V and the face opposite Vi, and the derivatives are 
taken in the direction of this same ray. See Figure 20.75. If we form the convex combina- 
tion 


5[F] 


f>]b'ibfSi[F] + bUlbfSj[F} + bjbfbj S k [F] + bplbfSdF] 

w+w a +w+w 

then 5[F] will match position and derivative values on the entire boundary of Tijki. 


C 1 , Discrete Interpolation in Tetrahedra 

For a C 1 , discrete interpolant, we assume that position and first-order derivative information 
is given at all four vertices of the tetrahedron Tijki. Since there are three (linearly indepen- 
dent) directional derivatives at each vertex, this amounts to a total of 16 data values. The 
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Figure 20.75: The data fora 16 parameter, Cl interpolant over a tetrahedron. 


method for describing an interpolant that will match these 16 pieces of data and which 
also has the property that all first-order derivatives across a face with common data will 
be continuous — is somewhat different from the previous interpolants we have described so 
far. Our description (and subsequent implementation) is based upon a two-step procedural 
discretization process. We use the transfinite interpolant of the previous section. In order 
to apply this transfinite interpolant, we need to define position and derivative values on the 
entire boundary of T ijk b First we assume that information is known on all the edges of the 
tetrahedra, and we describe how to extend it to the entire boundary. Second, we describe 
how to provide this transfinite edge data from only the discrete data at the vertices. If we 
know both position and derivative information on the edges, then we can use any C 1 trans- 
finite planar triangular interpolant to define position values on the interior points of the face 
triangles. For example, the side-vertex method itself could be used. Specifying position 
information on a face also implies some information about the derivatives on the interior 
of a triangle. Namely, all directional derivatives in a direction parallel to the face triangle 
are determined; so, in order to completely specify all derivatives, we need only provide a 
definition for the derivative perpendicular to the face. For this we use the C° version of 
the side-vertex interpolant which interpolates position data only and not derivatives, but we 
apply it to the edge data consisting of derivatives normal to a face. We now describe the 
second step of the discretization, which is how to compute edge information when only the 
point and derivative values are known at the four vertices. For position only on an edge, we 
simply use univariate cubic Hermite interpolation. This will also specify one directional 
derivative on the edge— namely which will vary as a quadratic polynomial. In order 

to get a C 1 join from one tetrahedron to the next, the other two directional derivatives must 
vary linearly along this edge. This is accomplished by specifying the gradient, VF, by the 
relationship 


= (l-t)VFi + tVFj 

' a p 

+ ^-(p)((l-f)VF I +tVF„e lj ) 

a e ij 


Zij 


V Fij (p) 


( 20 . 12 ) 
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where VFi = (F r (p,-) t F y (pi), F z {pi)) and t = |[j^~ p p j|| . This interpolation of the gradient 
is consistent with the value jf- already specified because (VF lJ (p), e^) = and it 
also has the property that for (n, e,j) = 0, 

(VFijip), n) = (1 - 0(VF,-, n) + t(VF Jr n), 


and so we have linear interpolation for any derivative in a direction perpendicular to 
eij. This completes the definition of the 16-parameter, C 1 , tetrahedral interpolant which 
is based upon the face-vertex interpolant. Examples and more discussion on this inter- 
polant can be found in [187]. The Clough-Tocher interpolant has been generalized to n- 
dimensional by Worsey and Farin [261]. Other C 1 , discrete interpolants for a tetrahedral 
domain are discussed in [2], [3], and [260], but each have some problem or drawback. The 
method of [2] is based upon the side-side, transfinite method of interpolation and appar- 
ently it has a problem with the linear independence of the discretized data. The method of 
[3] requires C 2 data for a C 1 interpolant and the method of [260] has a problem similar to 
its bivariate precursor [199] and [198]. This problem lies in the constraint that the center 
of the circumcircle of each triangle must lie interior to the triangular domain. 
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Abstract 

The interval volume is a generalization of the isosurface commonly associated 
with the marching cubes algorithm. Based upon samples at the locations of a 3D 
rectilinear grid, the me algorithm produces a triangular approximation to the 
surface defined by F(x,y,z) = c. The interval volume is defined by 
a < F(x, y,z)< p .We describe an algorithm for computing a tetrahedrization 
of a polyhedral approximation to the interval volume. 
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1. Introduction 


In this paper we describe an algorithm for computing a tetrahedrization of an interval 
volume. The interval volume is a generalization of the isosurface commonly associated 
with the marching cubes algorithm (see [4]). It is assumed that a trivariate function 
F(x,y,z) has been sampled at domain points lying on a 3D rectilinear grid. Given a 
threshold value, c , the me algorithm produces a triangulated approximation to the surface 
Sf( c ) = {(*>.y> z ) : F(x,y,z) = c} . Die interval volume is defined 

as If(a,/3) - {( x,y,z ): a < F(x,y,z ) < /?} . We describe an algorithm for computing a 
tetrahedrization of a polyhedral approximation to the interval volume. The uses and 
benefits of the interval volume have been delineated rather well in earlier papers on this 
subject (see Guo[3] and Fujishiro, Maeda and Sato [2]) and so we will not repeat them 
here. The algorithm of Guo [3] is based on the alpha shapes of Edelsbrunner and Mucke 
[1]. The algorithm of Fujishiro, Maeda and Sato [2] computes for each voxel the 
mtersecton of two convex polyhedra; namely I F (-co,/3) and I F (a,+cc) which are 

approximated by the polygon surface of the me algorithm [4] appropriately adjusted forht 
ambiguous cases (see Nielson and Hamann [6]). For the discussion here, the two surfaces 
which bound the interval volume, 5y(a)and S F (j3) are referred to as the a-surface and 
the /^-surface respectively. 



While it is impossible to anticipate all the application domains of algorithms for 
computing interval volumes, in this paper, we are primarily interested in applications 
where samples of the data function F(x,y,z) over a 3D rectilinear grid are known. This 
includes the typical data sets associated with MRI and CAT scans. In order to make the 
interval volume well defined, we make some inference about the variation of the data 
function, F(x, y, z ) , over the voxel domains defined by the vertices of the rectilinear grid. 
Our approach is very simple and leads to a two step algorithm: 
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Step 1. Each voxel is decomposed into tetrahedra and F is assumed to vary 
lineary over each tetrahedron. 

Step 2. For each tetrahedron of each voxel, the interval volume is computed and 
decomposed into tetrahedra. 

For the first step, there are several choices for decomposing a voxel into a collection of 
tetraheda. Since we require that the final collection of tetrahedra to be a proper 
tetrahedrization of the complete interval volume it is important that this first tetrahedral 
decomposition step also be a tetrahedrization. This implies, among other things, that the 
diagonals form the tetrahedization of one voxel to the next must match. We normally use 
the decomposition shown in Figure 2 which leads to five tetrahedra per voxel. The five- 
tetrahedron decomposition must be alternated from voxel to voxel with a rotated version 
of itself in order to maintain a proper tetrahedrization. 
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Figure 2. Steps of the Interval Volume Tetrahedrization Algorithm 


2. An Algorithm for Tetrahedrizing the Interval Volume 

The section is devoted to the details of the second step of our algorithm. At this point, 
we assume that we have a tetrahedrization of the domain of interest and the values of the 
data function, F, are known at the vertices of all the tetrahedra. We extend F to the 
interior of each tetrahedron by assuming that it varies linearly. With this assumption, the 
interval volume, I F {a,fS ) , is well defined and consists of a collection of polyhedron. Our 
algorithm computes the interval volume tetrahedization for each tetrahedron separately 
and so we need only describe the algorithm as it applies to a single tetrahedron with 
consideration to how the individual pieces will match up properly. The four vertices of 
each tetrahedron are classified according to the function value at this point. There are 
three possibilities. If the function value is less than or equal to a, then we call this a 
"white vertex" and mark it in our illustrations with an open white circle. If the function 
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value is strictly between a and j3, the we call this a "gray vertex" and mark it with a gray 
filled circle. If the value is greater than or equal to J3, then we call it a "black vertex" and 
mark it with a black filled circle. This classification leads to 15 distinct configurations 
which are all shown in Figure 3. All other cases can be rotated into one of these 15 
equivalence class representers. The labels we use for these configurations is based upon a 
triple index indicating the number of white, gray and black vertices present. This naming 
convention is shown in Figure 4. 





o < a 

a < © < p 

p < • 

Figure 3. The fifteen distinct configurations 
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# of o ^ y r ^ # of • 

# of • 

Figure 4. Method of labeling the configurations 

Note that all the polyhedra representing interval volumes in Figure 3 have planar faces 
and are convex. The planarity of the interior quadrilaterals lying on the a- or /^-surfaces is 

an immediate consequence of the use of linear interpolation to compute the vertices along 
the edges. 6 

We next describe how to decompose each of the interval volumes of Figure 3 into 
tetrahedra. We first note that there are four distinct types of polyhedra, namely: 

Tetrahedra: Configurations: 013, 040, 310. 

Prism shaped polyhedra: Configurations: 031, 022, ,103, 130, 220, 301 

Each of these polyhedra has six vertices, two opposing triangular faces and three 
planar quadrilateral faces. 

Crystal shaped polyhedra: Configurations: 1 12, 121, 21 1 

Each of these polyhedra has eight vertices, two opposing triangular faces, two 
planar quadrilateral faces and two planar pentagonal faces. 

Cube shaped polyhedra: Configuration: 202 

The polyhedra of this configuration has eight vertices and six planar quadrilateral 
faces. 

We now describe how each of these distinct types of polyhedra can be tetrahedrized so 
as to lead to a global tetrahedrization of the interval volume. Since they are all convex 
polyhedra, there is no problem as to whether or not they can be individually tetrahedrized. 
(See Nielson [5].) The problem becomes interesting in trying to maintain a global 
tetrahedrization. Unless special provisions are made, it is possible that the tetrahedrization 
of a portion of the interval volume in one tetrahedra could not join with the 
tetrahedrization of another portion of the interval volume in a separate tetraheda in manner 
that would lead to an overall tetrahedrization of the interval volume. This could happen, 
for example, if these two portions shared a common quadrilateral face and in one instance 
one diagonal is chosen and in another instance the other diagonal is chosen. See Figure 5. 
This eliminates the possibility of having a well-defined piecewise linear function defined 
over this type of decomposition. In order to circumvent this problem and have a consistent 
choice of edges on any of the planar faces, we use the “index connection rule.” The 
invocation of this rule assumes that we have a total ordering on the vertices so that for any 
two vertices one is "<" or "smaller" than the other. There are, of course, many ways to 
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accomplish this linear ordering of the vertices. The details of how we do it for the 
examples presented here is given in the next section. 




Figure 5. The "3D cracking problem" on the right. 

Index Connection Rule: 

Quadrilateral: The edge consisting of the unique smallest vertex joined to its 
opposing vertex will be contained in the tetrahedrization. 

Pentagon: Each pentagonal face has a unique vertex which is gray; that is, 

the function value at this vertex is between a and J3. Also, it is 
always the case on pentagonal faces that there are two adjoining 
a-vertices and two adjoining /^-vertices. The gray vertex and its 
non neighbor a-vertex will form an edge as well as the gray 
vertex and its non neighbor /7-vertex. See Figure 7. 



Figure 6. Index connection rule applied to quadrilateral. 
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Figure 7. Index connection rule applied to pentagonal face 


We now set out to prove how each of the different types of polyhedra (prism, crystal 
and cube) can be tetrahedrized and in a manner which is consistent with the edges required 
by the “index connection rule.” 

Prism shaped polyhedra: Configurations 031, 022, 103, 130, 220, 301. On any prism 
shaped polyhedron, there are eight different possible edge connections for the three 
quadrilateral faces. Six of them can be “realized” by a tetrahedrization of the polyhedron 
while two can not. See Figure 8. The edge connections shown in the upper left and lower 
right are the ones that can not be realized, but these particular vertex connections would 
never be specified by the "index connection rule" or otherwise we would have by 
transitivity a contradiction to the linear ordering of the vertices. 
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Figure 8. The eight possible edge configurations for the quadrilateral faces of the prism 

shaped polyhedron. 

Crystal Shaped Polyhedra: Configurations: 121, 112,211. For this type of polyhedra, 
we can first decompose it into two pyramids and one tetrahedron and this decomposition 
leads to edges on the boundary of the polyhedron which are consisitent with the edges 
specified by the “index connection rule” applied to pentagonal faces. This is shown for 
the three configurations, which lead to this case, in Figure 9. The “index connection rule” 
is now applied to the remaining quadrilateral faces of the pyramids and these are split into 
two tetrahedra each leading to a total of five tetrahedra decomposing a crystal shaped 
polyhedron. 
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Figure 9. Decomposition of Crystal shaped polyhedra into a tetrahedron and two 

pyramids. 

Cube Shaped Polyhedra: As with the other cases, we must show that the polyhedron 
can be tetrahedrized such that the edges imposed by the “index connection rule” are 
realized. Not all possible edge configuratons on the faces of this hexahedron can be 
realized by a tetrahedrization. See Nielson [5] for examples that are not realizable. But 
the edges imposed by the “index connection rule” are all fully realizable as we now will 
show. First, let us identify the smallest vertex of this polyhedron and also the second 
smallest. In the figures, the smallest is marked with a "1 " and the next smallest is marked 
with a "2". There are three cases to consider. The first is when these two special vertices 
are located diagonally opposite from each other. Since everyone of the eight faces 
contains exactly one of these two special vertices, the “index connection rule” applies and 
the edges on all eight faces are completely specified. Fortunately, this particular 
configuration of edges is fully realizable by a tetrahedrization of the cubed shaped 
polyhedra. The tetrahedrization consistent with these edges on the faces is shown in 
Figure 10. 


1 

Figure 10. The first type of case for the cube shaped polyhedron. 

The next case is when the smallest vertex and the next smallest vertex are diagonally 
located on the same face. All faces except the face opposite the one containing these two 
special vertices has at least one of these special vertices and so the edges are determined 
for these faces. This leaves ony two cases to consider. These two cases are illustrated in 
Figure 11. In the first case, shown on the left, all opposing faces have diagonal edges 
which are switched. This particular edge configuration can be realized with a 
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tetrahedrization consisting of five tetrahedra as shown. The other case shown on the right 
is realized with the tetrahedrization shown consisting of six tetrahedra. 



The third case is when the smallest vertex and the next smallest are on a common 
edge. In this case there are four subcases as the “index connection rule” leaves the edges 
on two adjoining faces undetermined. These four subcases and the realizing 
tetrahedrizations are shown in Figure 12. In all subcases, we have a tetrahedrization 
consisting of six tetrahedra. 
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3. Implementation Considerations and Examples 

Our implementation produces a tetrahedral grid rather than simply a list of the four 
vertices of each tetrahedron. The efficiency and benefits of a tetrahedral grid over this 
simpler data structure are evident for most applications. For compactness in the 
following discussion, rather than the cumbersome "tetrahedral grid representing a 
tetrahedrization", we often use the simpler "tetrahedrization." Our definition of a 
tetrahedrization starts with the collection of vertices P t = (x h y h zi), i = 1 which 
we assume are not collectively coplanar. We denote this collection of points by P . A 
tetrahedrization consists of a list of 4-tuples which we denote by /, . Each 4-tuple, 

ijkt e I, denotes a single tetrahedron with the four vertices Pj,Pj,P k ,P e . The following 
conditions must hold: 

i) No tetrahedron T ijkh ijkl g I, is degenerate. That is, if ijkl e I t then 
Pi, PjP k and P ( are not coplanar. 

ii) The interior of any two triangles do not intersect. That is if ijkl e /, and afiyS g I t 
then \nt(Tjj k p ) n Int(T^) = <j> . 

iii) The boundary of two tetrahedra can only intersect at a common triangular face. 

iv) The domain is the union of all tetraheda, D = U T ijkl ■ 

ijkt el : 

We should point out that condition iii) must hold in the strictest sense and so 
tetrahedra joining as shown in right side of Figure 5 are not allowed. The reason for this 
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condition (and all the others) is that we eventually wish to be able to define C® functions 
in a piecewise linear manner over the domain consisting of the union of all tetrahedra. 

The tetrahedral grid data structure for representing tetrahedrizations is illustrated by 
the example of Figure 13 where a tetrahedrization of the cube into 5 tetrahedra is shown. 


Y 



Figure 13. An example which defines the tetrahedral grid data structure. 


As we noted earlier, the use of the "index connection rule" requires that we have a 
linear ordering of the vertices. The ordering we use is based upon the lexicographic 
ordering of a unique naming rule. The name of a vertex is specified by 

vertex_name = i, j, k, X|Y|Z|XY|XZ|YZ|YX|ZX|ZY|0,a|p (1) 

This is further illustrated in Figure 14. As the algorithm processes each tetrahedron (five 
per voxel) this naming scheme is used to record the tetrahedra of the interval volume for 
each tetrahedron. This makes it easy to invoke the "index connection rule" locally. Later 
we post process this tetrahedral grid into the more conventional format illustrated in 
Figure 13 by sequentially replacing the name (pointer) as describe in equation (1) with a 
simple integer pointer (or one-dimensional array index). It is also interesting to note that 
this particular ordering eliminates the possibility of the first and second types of cases for 
the cube shaped polyhedra and so only the third type of case (Figure 12) must be 
considered. 
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i, j, k X i+l,j, k 

Figure 14. Vertex naming scheme. 


Example 1. 

This first example utilizes CAT data on a 64 x 64 x 68 grid. The lower value for the 
interval volume is a = 50.0 and the upper value is /? = 150. There are 362,181 

tetrahedra and 115,956 vertices. In Figure 15 we show the interval volume. The a- 
surface is colored white and the /2-surface is colored red. We have cut a section of the 
data away so as to show the interior of the volume. Triangular faces of the interval 
volume which do not lie on either the ar-surface or the /2-surface are colored blue. 



Figure 15. Display of interval volume for Example 1. 

Example 2. 

This next example is based upon some MRI data. The grid size is 64 x 64 x 58 and 
we have chosen the a and (3 values in order to try and get a volume representation of the 
"skin". These are taken to be 15.0 and 40.0 respectively. These values lead to 454,916 
tetrahedra and 147,71 1 vertices. A rendering of the interval for this example is shown in 
Figure 16. The left image is the interval volume and we see nothing more than the a- 
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surface which has been colored white. In the right image we show a cut-away with the 
triangular faces not lying on a either the a-surface or the ^-surface colored in blue. 



Example 3. 

The next example is based upon the data function 


F(x,y,z) = 12 y + 36(x-z) 2 -5 


and the values: 


a = -0.543, p = 0.0 , 

Domain: {(x,y,z): 0 < x < 1, 0<y<l, 0<z<l} 

Sample Grid: i = 0.....14; j = k = 0,...,14}. 

The final tetrahedrization has 4,034 tetrahedra and 1,480 vertices. A picture of the interval 
volume is shown in Figure 17. 
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Example 4. 

This example utilizes the data function 


F(x,y,z) = 4(3 y - l) 2 + 18(x - z) 2 - 18(x + z - 1) 2 + 1 


and the values: 


a = - 1 . 1 1,/3= 1.1 1 


Domain: {(x,y,z): 0<x<l, 0<y<l, 0<z<l} 

Sample Grid: / = 0 ,...,14; j = 0.....14; k = 0.....14} . 

The interval volume consists of 5,782 tetrahedra and 2,092 vertices. A rendering of the 
interval volume is shown in Figure 18. 
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Figure 18. Interval Volume of Example 4. 


Example 5. 

This example is based upon the spherical data function 


F(x, y, z) = (x - 0.5) 2 + (x - 0.5) 2 + (x - 0.5) 2 


and the values: 


a = 035,0 = 037 

Domain: {(x,y,z): 0<x<l, 0 < jv < 1, 0<z<l} 

Sample Grid: {(^,^,^>, i = 0,...,14; j = 0,...,14; k = 0,...,14). 

The tetrahedrization has 8,500 tetrahedra and 3,224 vertices. A rendering of the interval 
volume is shown in Figure 19. 
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Figure 19. Interval Volume of Example 5. 


Example 6. 

The data function for this example is 


F(x, y, z ) = a(x - 0.5) 2 + a{y - 0.5) 2 + {a - 1 )(z - 0.5) 2 , a = — 

16 

and the values of the other parameters are: 

a = -0.02, p = 0.02 


Domain: {(x,y,z): 0<x<l, 0 < jv < 1, 0<z<l} 

Sample Grid: i = 0,...,14; j = 0,...,14; k = 0,...,14}. 
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Figure 20. Interval Volume of Example 6. Shaded and wire frame renderings 


In this example, we have also included a wire frame renderings so that the 
tetrahedrization is apparent. 

It is interesting to note the frequency of the various configurations. In Figure 21, we 
show the statistics for the frequency of occurrence of these configurations for the six 
examples we have discussed here. 
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4. Remarks 


1 . While our primary attention for application data sets has been rectilinear grids, for the 
most part, our algorithm works for any tetrahedrization of a domain. This includes 
tetrahedrized 3D curvilinear grids and scattered volumetric data sets. 

2. Some readers may wonder why we did not use the Delaunay tetrahedrization of the 
convex polyhedra of Step 2 since a Delaunay tetrahedrization would lead to a two- 
dimensional Delaunay choice on the common faces and therefore gaurantee a matching 
across common faces. The problem here is to properly and efficiently deal with the 
neutral and near neutral cases which are far from uncommon in the applications that 
we have discussed here. 
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Abstract 

Some new piecewise constant wavelets defined over nested 
triangulated domains are presented and applied to the problem of 
multiresolution analysis of flow over a spherical domain. These new 
nearly orthogonal wavelets have advantages over the weaker biorthogonal 
wavelets. In the planar case of uniform areas, the wavelets presented here 
converge to one of two fully orthogonal Haar wavelets which are proven to 
be the only wavelets of this type. 
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1 Introduction 


In this paper we present some new Haar (piecewise constant) wavelets defined over 
nested triangular domains. Our target application is the multiresolution approximation 
and analysis of fluid flow over a spherical domain and we will show the application of 
these new wavelets in this context in Section 3. Even though our primary application is a 
spherical domain and for simplicity, our figures will be for the case of planar triangles it 
should be kept in mind that the development and application is much more general. The 
techniques developed here apply to any triangulated manifold or surface. 


The concept of the orthogonality of the wavelet basis functions is important. This 
relates to the ease of computation of best approximations. Fully orthogonal (or simply 
orthogonal) wavelets are to be preferred over the weaker biorthogonal in this regard. 
ITiere do not exist orthogonal wavelets over spherical nested triangular domains (only 
biorthogonal) and this paper does not change that. However, the new wavelets presented 
here are "nearly orthogonal" in the sense that when the domain becomes planar (as 
spherical grids approach) the wavelets become fully orthogonal. Exactly two, planar 
a me invariant, orthogonal wavelets result as limits of our nearly orthogonal wavelets 
and we eventually prove that these are the only planar wavelets of this type. 

Often these days, papers on wavelets and multiresolution approximation and analysis 
will give some general basic motivation and overview. We will dispense with that here 
and refer the reader to anyone of the many excellent introductions currently available. 
See 11 J for example. Rather, we will focus our introductory material on the case of Haar 
wavelets over nested triangular grids. First, some comments about the types of domains 
we are covering. We start with an initial triangulation of a domain of interest This 
mitial or base triangulation may be very simple. In the case of a planar application it may 
consist of a single triangle. In the case of a sphere it may consist of the eight spherical 
triangles with vertices on an octahedron or four spherical triangles with vertices on a 
tetrahedron. Each of the triangles of the base triangulation is subdivided into four 
subtnangles in a manner shown in Figure 1 by adding a vertex along each edge Often 
this vertex will be the midpoint (geodesic bisector), but this need not be the case A 
generic triangle of the base triangulation is denoted by Tand we concentrate our attention 
from this level forward. All subtriangles are denoted by T n , where n is an N-tuple made 

from the symbols o, i, j or k. The triangle labeling scheme that we use is illustrated in 
Figure 1 . This labeling scheme allows us to dispense with the normal data structure used 
to represent a triangular as the name itself allows for the determination of all neighbors. 

The parent of T o , T ni , T nj and T nk is T n . At level N there are 4^ subtriangles 
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Figure 1. The Triangle Labeling Scheme for Subtriangles 


We start with a function which is defined in a piecewise manner over a certain level 
of the nested triangular domains and we wish to find an approximation defined in a 
piecewise manner over the coarser, parent grid one level up. We usually think of 
applications starting with a rather fine resolution so that simple, piecewise constant 
approximation is sufficient. We want to compute this lower resolution approximation in 
an efficient manner and in a manner that allows the reconstruction of the original 
piecewise defined function. This is accomplished through a change of basis which we 
now illustrate by example. We can represent an arbitrary piecewise constant function 
over the triangular subdivision at level 2 (upper right of Figure 1) by 

F (p) = hkhk ( P ) + hofacAP) + hihAP) + KjAAp) + Aj</>ij(p) + A o 0io(p) + h<t>ii(p) 

+ ^kifikiiP) + A-oifioAP) + ^oo^oo(P) + ^okfiojiP) + Ak^ik(P) 

+ Ajitji(p) + AjotjoiP) + Aj k <Pjk(P) (1 

+ A.jj<f>jj(p). 


We have arranged the terms of F(p ) so as to correspond with the layout of Figure 1. 
The value A ab represents the value of F(p) over T ab . The first step of the so called 
wavelet decomposition process yields 
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F (p) = rkkVjjip) + hh(p) + YkjVkjip) + y OJ Voj(p) + YijYij(p) + ^Mp) + YuVuip) 
+ Yki'PkiiP) + YoiVoi(p) + *vt 0 (P) + YokVokiP) + YikYik(P) 

+ Y jiVji(p) +Aj</>j(p) +Y jk Vjk(p) 

+ YjjVjj(p) 


where the coarser approximation is represented by the terms 


h$k(P). , j. , n s + Wi(P) 
+ AAO) 

+ A j<pj(p) 


and the so called detail or wavelet is represented by the 


remaining 12, if/ basis functions. We should be quick to point out that F(p)has not 
changed from equation (1) to equation (2), only a change of basis has taken place. Still F 
has the property that F{p) = \ ab for p € T ab . We can carry this basic idea one step 


further and replace the coarser approximation, + X 0 <j> 0 (p) + , which is 

+ kjfijip) 

piecewise defined over the parent triangles with an even yet coarser approximation which 
is, in this second and last step, piecewise over the single triangle parent of 7), Tj,Tk,T 0 . 

This leads to the representation 


F (P) = Ykk ¥kk(P) + YkFk(P) + YkjYkj(P) + YoiFoi(P) + YijVij(p) + YiVi(p) + YiiFn(p) 
+ YkiVki(P) +Yoi'Foi(p) + X(f>{p) +Yok'Fok(P) + Yik'Pik(P) 

+ Yji*Fji(P ) +Yj¥j(P) +YjkVjk(P) 

+ YjjVjj(p). 


From equation (2) to equation (3) the 12 wavelet terms from equation (1) remain and 3 
additional wavelet terms have been added along with the coarser (coarsest) approximation 
X<f>{p) ■ A general step in this decomposition is accomplished by successively processing 
groups of four terms which are associated with a particular triangle. The terms 

^nkfinkiP) ^-nofino^P) ^7ii t fini(P') 

+ ^nj<f>nj(P) 


are replaced by the terms 


using the relationship 


Ynk'Fnk(P) + Ki<t> n {p) + Yni<Fni(P) 
+ YnjFnjiP) 
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The above relationship (4) is called the refinement relationship within the context of 
wavelets and embodies the information for a change of basis functions. See Figure 2. 
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Figure 2. Notation and an Illustrated Representation of Wavelet Decomposition 


Actually in the computation of the coefficients of the wavelet decomposition leading 
from (1) to (3) we need the inverse of the refinement relationship; namely, 




r Ko 

Snoi 

Snoj 

Snok 

r UP)" 

<Pni(P) 


Ki 

Snii 

Snij 

Snik 

¥ni(P) 

<Pnj(p) 


Kj 

Snji 

Snjj 

Snjk 

¥nj(P) 

\UiPP 


K h nk 

Snki 

Snkj 

Snkk) 

Wnk(P)J 


because 


5 




( 6 ) 


(^77 Yni Y nj Y nk) (^770 ^ni ^ nj ^ nk ) 


^no 

Snoi 

Snoj 

Snok 

Ki 

Snii 

Snij 

Snik 

Kj 

Snji 

Snjj 

Srijk 

^nk 

Snki 

Snkj 

SnkkJ 


This inverse relationship is called the decomposition relationship. The representation of 
equation (3) is called the wavelet decomposition and from this representation, we can 
easily reconstruct any of the piecewise defined approximations at any level and this can 
be done either locally or globally. Also, from this decomposition, we can easily 
reconstruct the original function locally by simply adding up the terms whose support 
intersects the local region of interest. Other oracles can be invoked for different purposes 
and applications. An oracle is a means by which a certain subset of the basis functions of 
the wavelet decomposition are selected to form a certain beneficial approximation of the 
original function. Compression of the original function is a very popular application. 
Here a certain subset (say those with the largest coefficients) of the terms of the wavelet 
decomposition are maintained. If the coefficients of the discarded terms are nonzero, 
then there will be some loss as a result of the compression. With regard to this, the 
concept of orthogonality is important. Let us explain why by noting a very basic property 
of best approximations and orthogonal basis functions. Suppose that we wish to solve the 
problem of 


min 

{^1 , CZ 2 > ' * ’ j } 


N 

F-X a iSi 


/'=! 


where ||F|| 2 = (F,F). If the basis functions are orthonormal, i. e. = 

then the solution of (6) is easily computed as 

N 

7=1 

Furthermore this implies that whenever you have a representation of a function of the 
form 


1 1 * = j 

[ 0 , i * J 


K 

F = H a iSi 

7 = 1 

where the basis functions are orthonormal, then a, must be (F,g, ) and if we wish to 
compute a best approximation to F from among a subset of these basis functions, we need 
only keep those of interest. That is, if F is given by (7) then the solution of 
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min 


a h i eK s 


Tu a iSi 

i*K. 


is 


G = X a i§i ■ 

isK, 

This only works if the basis functions are orthogonal. 

Biorthogonal wavelets have the property that the collectively the y/’s are orthogonal 
to the (jf s , but individually, the y/' s are not orthogonal to each other and so by truncating 
a wavelet representation, we don't necessarily get the best approximation. As long as we 
do the approximation in "chunks" it is OK. That is, as long as we involve all the wavelet 
basis functions from a particular level the best approximation property holds, but if we 
use an oracle which does not keep this in mind, then we are not making the best use of the 
storage for the coefficients. For example, an approximation based upon a certain 
percentage of the largest coefficients of the biorthogonal wavelet decomposition will not 
necessarily give the best compression in the least squares sense. 

There is a second concept that is really quite important when it comes to the practical 
application of these techniques. This has to do with how the actual data is mapped to the 
names and variables used to describe and program these methods. For example suppose 

that we start with a triangular array of trixel data Fy, i,j = l,---,2 ,/ + j <2 N and we 

plan to apply the wavelets described here to obtain some compression or lower resolution 
approximation. In order to apply these techniques, we need to make some association of 
the data to the subtriangle domains illustrated in Figure 1 . There are six possible ways of 
making this association with each being an affine map of the other. But there is no 
particular reason why one would be preferred over the other. We certainly do not want 
the results of our wavelet decompositon and expansions to depend upon which choice is 
made. Therefore, it is important that are methods give the same results in any case and be 
invariant to these affine transformations or equivalently be invariant under the labelings 
used. Some authors have called this property "symmetry." 

2. Examples of Haar Wavelets over Triangular Domains 

Within the context that we have established here, a particular type of Haar wavelet is 
completely determined by the refinement equations of (4). We begin this section by 
presenting two previously discussed Haar wavelets. The first is mentioned in [3] where it 
is attributed to Mitrea and Girardi & Sweldens. The refinement equation for these 
wavelets is 
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where a n represents the area of triangle T n and a ni ,a n j,a nk ,a no represent the areas of 
the subtriangles T nj , T n j , T nk , T no respectively. These Haar wavelets are orthogonal in 
that y/ ni , y/ nj , y/ nk are all orthogonal to <f> n and also y/ ni ±y/ nj , y/ ni Ly nk , y/ nj Ly / nk . Also 

JVm| = J| ¥nj \ ~ \\¥nk\~ a - Unfortunately, these wavelets are not affine invariant 
T, r T, 

which means that the results of any application of these wavelets would depend upon the 
happenstance of how the vertices were labeled. 

Another type of wavelet (see [3]) is defined by the refinement equation 
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These wavelets are affine invariant, but this comes at the cost of giving up the 
orthogonality. These wavelets are only biorthogonal. Also, as above, 

— ^\¥nj | = JW| = «. 


We now present our two new types of nearly orthogonal wavelets. The first has the 
refinement equation 
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and the second has the refinement equation 
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Both of these new wavelets are biorthogonal and affine invariant but not orthogonal. But 
these new wavelets have an advantage over other wavelets of this type. As the areas of 
the subtriangles become more nearly equal, these wavelets tend to fully orthogonal 
wavelets. This is not true of the wavelets given by equation (7) for example. If we let 
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and so we have the innerproducts: j \¥ni¥nj = \vniVnk ~ j ¥nk¥nj = In the case 
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of the nearly orthogonal wavelets, we have, for uniform areas, the refinement equations 
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and it can easily be verified that in both cases <j> n 1 y/ nh <f> n _L y/ n j, <j> n _L \(/ nk and 

Wni -L Wnj » Vni -L nk > V nj -L Vnk • So, in both cases, as the areas of the subtriangles 

become uniform the wavelets are become mutually orthogonal. It is interesting that these 
two wavelets are the only affine invariant, fully orthogonal Haar wavelets. 

Theorem: In the case where the areas of the subtriangles are all equal, there are only two 
types of fully orthogonal, affine invariant Haar wavelets; namely those defined by the 
refinement equations given in equations (11) and (12). 

Proof: If we impose the affine invariance conditions then we must start with refinement 
and decomposition relationships of the general form 
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Now if we further impose the conditions of full orthogonality, we have the requirement 


(gg-g o )Oge + 5g o ) = 0, 


which give rise to two solutions. The first solution (g e = g 0 ) leads to the wavelets of 

equation (11) and the second solution (g e = -|g 0 ) leads to the wavelets of equation 

(12) and so these are the only, orthogonal, affine invariant Haar wavelets over subdivided 
triangular grids. 

3. Application Examples 


We have applied our new techniques on a variety of data sets, but some of the most 
interesting are for the case of a spherical domain. We present some of those results in 
this section. First we mention a few things about the triangular grids for a spherical 
domain. The construction of the grid starts with a base triangulation of the sphere. Often 
this is based upon the vertices of a tetrahedron (4 triangles) , octahedron (8 triangles) or 
icosahedron (20 triangles). Each triangle of the base triangulation is subdivided into four 
triangles by bisecting the geodesic edges at their midpoints. This process is continued. 
Figure 3 illustrates this process with the base triangulation being a spherical tetrahedron. 
Even though the triangles of Figure 3 are drawn with straight edges, the actual triangles 
used for the wavelets are spherical triangles with geodesic arcs as edges. For the sphere, 
it is impossible to arrange for the areas of the subtriangles to be equal while in fact the 
areas of these subtriangles can be quite different. For example, in the first level of the 
case of a tetrahedron as the base triangulation vertices, the area of the center subtriangle is 
four times the area of the the other three subtriangles! As we go deeper into the 
subdivision, this ratio diminishes to one. It is this aspect of spherical triangles that is 
exploited by the nearly orthogonal wavelets we have presented here in equations (8) and 
(9). Because of this significant difference in the areas for the tetrahedron at lower levels, 
we usually use a base triangulation that utilizes the vertices of an octahedron and this 
includes the examples presented here. 

We have implemented the wavelets of equation (7) and the new wavelets of (8) and 
(9) and applied them to a variety of spherical data sets. For spherical images, the results 
of Table 1 are typical. Here we took a spherical image at level 7 ( 32,768 trixels) and 
computed the wavelet decomposition and then compared the RMS errors of the 
approximations using 10%, 50% and 90% of the terms with the largest normalized 
coefficients. Most generally the nearly orthogonal wavelets out performed the wavelets 
of (7), but there does not seem to be a clear choice between the two different nearly 
orthogonal wavelets. 
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Level 3 


Level 6 


Figure 3. Spherical Triangular Grid 


More interestingly, we have applied our new wavelets to data consisting of a vector 
field defined over the sphere. In the example of Figure 4, we took a known vector field 
function and evaluated it to obtain the data. We applied the wavelet of equation (8). In 
the left column we show three views of the reconstruction of the wavelet decomposition 
using only 1% of the coefficients. This means we started with a vector valued function 
which had 2048 coefficients (one for each triangle). The first time we apply the 
decomposition equations we obtain a function which is piecewise constant over 512 
triangles and plus the 1536 detail (or wavelet) functions. No information is lost, just a 
change of basis. The next step yields a lower resolution approximation consisting of 128 
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terms and 384 detail at this level plus the 1536 detail from the first decomposition step. 
This is continued two more times. This final expansion (of 2048 terms) is scanned and 
the 20 (1%) functions with the largest coefficients are summed to obtain the partial 
reconstruction. In the right column we show the reconstruction based upon the largest 
409 (20%) coefficients. It is indistinguishable from the original function. It is interesting 
to note that the main features of the flow are maintained with an approximation which 
uses only 20 out of 2048 terms and the flow an be compressed to 409 terms (20%) with 
no visual loss. The graphs of Figure 4 and Figure 5 are computed using the tangent curve 
advection and topological methods described in [2]. 
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View 3,1% View 3, 20% 

Figure 4. Partial (nearly orthogonal) spherical wavelet reconstruction of a flow field 

defined over a spherical domain. 
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In this next example shown in Figure 5, we applied the wavelets of equation (8) to 
some "real" data. This is data we obtained from Roger Crawfis and Nelson Max of 
Lawrence Livermore National Laboratory. Actually it is simulated data resulting from a 
global weather model. This data represents one slice in elevation and one time step. The 
first column shows three different views using a topological graph of the original data. In 
the second column the reconstruction based upon the 3% largest coefficients is shown. 
This low resolution approximation allows us to get an overview of the flow without the 
clutter and detail. The use of these types of models for analysis and visualization of 
scientific data is potentially very useful and a very exciting area of research. 
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View 3, 100% View 3, 3% 

Figure 5. Spherical wavelet reconstruction of wind data from global weather model 
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Abstract 

An algorithm for computing a triangulated surface which separates a 
collection of data points that have been segmented into a number of different 
classes is presented. The problem generalizes the concept of an isosurface which 
separates data points that have been segmented into only two classes: those for 
which data function values are above the threshold and those which are below the 
threshold value. The algorithm is very simple, easy to implement and applies 
without limit to the number of classes. 
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1. Introduction and Algorithm 

In this paper we describe an algorithm for computing a triangulated surface which 
separates regions of different types. We assume that we have a collection of data points 
V n i = 1 and that each of these data points has been classified into one of several 
possible classes c,,i = . This includes, for example, medical scanning device data 

that has been postprocessed by some segmentation procedure into different tissues or 
organ classes or physical simulation data that has been classified by material properties 
such solid, liquid or gas. For our algorithm, we assume that the data points V t - are the 
vertices of a tetrahedrization of the domain of interest. Two important application areas 
are where the data points lie on a 3D rectilinear grid or 3D curvilinear grid. In either case 
we preprocess the data by subdividing each voxel or hexahedron (curvilinear grid cell) 
into tetrahedra and proceed with our algorithm. See Figure 1 and Nielson [4], Note that in 
the 6 tetrahedra split, each cube is split exactly as shown. In the 5 tetrahedra split, the 
mirror image of the split alternates in each direction with the one shown. 





Figure 1 . Decomposing voxel data into tetrahedra data. 


Our goal is to produce a triangulated surface which separates the components 
(connected subsets) of the regions, R n i = each containing the data of type 

Cj, i = 1,..., M . This surface can be viewed as a generalization of the isosurface often 
associated with the marching cube algorithm (see [3]). In the context of the me algorithm 
the discrete vertices lying on a 3D rectilinear grid are classified into only two possible 
classes: either the value of the data function, 5 , is above the threshold of the isosurface or 
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below this threshold. The isosurface then separates these two classes of data points into 
two regions. In the more general situation where there are several possible classes for data 
points, the separating surface is defined as S = n . In the spirit of the me 

algorithm, our algorithm processes each tetrahedron separately. Let be the 

vertices of an arbitrary tetrahedron. If two vertices, say V { and Py are classified 
differently, we make reference to a point my along the edge joining V i and Vj . This 
point is where this edge intersects the surface separating the vertices Vj and Vj . This 

separating point can be anywhere on this edge and in some default situations a reasonable 
choice would be the midpoint. We mention some other considerations for selecting this 
point later. If three points of a face, ,Vj and V ) c , are classified differently, we must 

make reference to a point lying somewhere interior to this face. Again, for the 

topological aspects of our algorithm, it is not important where exactly this point lies on the 
face, but some practical considerations which we discuss later lead to reasonable choices 
for this point. And finally if all four points are classified differently, we need to reference 
a point m t lying interior to the tetrahedron. This notation is further illustrated in Figure 2. 



Figure 2. Notation used for vertices, mid-edge, mid-face and mid-tetrahedron points. 

The strength of our tetrahedral-based algorithm is its simplicity and subsequent ease of 
implementation. There are only five cases to be considered: (0) all vertices are classified 
as one type (trivial case; no separating surface intersects the tetrahedron), (1) three 
vertices are of one class and one other vertex is of another class, (2) two vertices are of 
one class and two vertices are of another class, (3) two vertices of one class and the other 


3 



two vertices are of second and third classes, and (4) each vertex is of a different class. 
Because any configuration in one of these five cases can be rotated into a standard 
configuration, standardized algorithms can be used assuming that (local) vertices are 
labeled V h Vj,V k and Vj. 

The face of a tetrahedron having vertices of more than one type must be split. This can 
be seen in Figures 3, 4, 5 and 6 for the four nontrivial cases. When the vertices on a 
particular face of the tetrahedron are of only two types, the face is split along the line 
segment joining the mid-edge points on that face, say the points and mjj . This occurs 

in cases (l)-(3). When the vertices on a face are all three of a different type, the face must 
be split not only at the mid-edge points, but also at the mid-face point interior to this 

face. The face is then divided by the line segments joining the mid-face point to the mid- 
edge points. When four different types are present then we must involve the mid- 
tetrahedron point m t . The separating surface is to be represented as a union of triangles, 
so quadrilaterals that naturally occur in our algorithm must be triangulated by including 
one diagonal or the other. We adopt the convention that we will impose those diagonals 
that are consistent with a certain tetrahedrization of the four hexahedra that occur in case 
(4). See Figure 6. Since each hexahedron has a vertex of the cube as one vertex, we adopt 
the triangulation of the faces by diagonals from the tetrahedron vertices to the mid-face 
points, and the mid-edge points to the mid-tetrahedron point. We should point out that 
unless certain restrictions are put on the mid-edge, mid-face and mid-tetrahedron points, 
those quadrilaterals may not be planar. This causes no particular problem, although we 
note that the separating surface will be slightly different if different choices were made 
when triangulating those quadrilaterals. 

In case (3), there is a dilemma as to whether the exterior quadrilateral faces 
(V h V j,mji,mji and Vj , V j , mj k , m ik in Figure 5) would be divided from V i to mjj and V i 

to mjk , or from Vj to and V j to when tetrahedrizing the subvolumes. We have 

adopted the strategy that the order- of the vertices in the description of the tetrahedrized 
volume will determine that; we choose Vj or Vj according to which has priority in the 

input list. Because we maintain order when sorting the unique classes of vertices for a 
particular tetrahedron, symbolically the first vertex is Vj. Hence we diagonalize the 
interior separating quadrilaterals using the line segments through and , and mj^j 

and mjk . 

In case (2), we again wish to be consistent with some tetrahedrization of the volumes, 
two triangular prisms in this case. Thus, we follow the rule of connecting the priority 
vertex to the opposite mid-edges for each prism, i.e., Vj to mjk and mjj , and to am,-/ 

and mjj . After this is done for each prism, it is seen that the diagonal on the separating 
quadrilateral is arbitrary, and we choose the mjk to segment. 
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For completeness we list the triangles comprising the separating surface in each case. 


Case (1): 

Triangle: m ih mj h m kl 


Vi 



Figure 3. Case (1): Three vertices of one class and one vertex of another class. 
Case (2): 

Triangles: m ik , m jk , m n ; m a , m jk , ntji 



Figure 4. Case (2): Two vertices of one class and two vertices of another class. 
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Case (3): 

. m kl > m ikl > Mjkl , m jk > m ik > m jkl > m ik > m jkl > m ikl > 

Triangles: 

m jl > m il > m jkl » w ;7 > m ikl > 171 jkl 



Figure 5. Case (3): Two vertices of one class and two other vertices each of another class. 


Case (4): 

m ij > m ijk > m t ; "ty . m ijl » m t ; m jk » m ijk » , mjkl , m t ; 
Triangles: m ik , m ijk , m t ; m ik , , m, ; , m j kl ,m t \m jh m tjl , m t ; 

m ;7 » m ijl » m t ; W// , m M , m t ; , w //t/ , m t ; m kl , m jkl , m t 



Figure 6. Case (4): Each vertex is a different class. 
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For the general description of our algorithm, we have kept the location of the mid- 
edge, mid-face and mid-tetrahedron points arbitrary. It is easy to present this way and also 
this allows for maximum flexibility. In some applications where there is no additional 
information on which to base any bias or adjustment, one just as well select these point to 
be the actual geometric midpoints. That is, 


Vi + V j 

171 :: = 

9 2 


m ijk = 


m ij + m ik + m jk 


m, = 


m ijk + m jkl + m M+ m ijl 

4 


( 1 ) 


In some other applications where there is additional information some weights may be 
used to compute these values. For example, if data points are classified (or segmented) on 
an interval of values of some data function, then it might be useful to weight accordingly 
the computation of the mid-edge value. Assume that an arbitrary point V = (x, y, z ) is 
classified by the rules: 


and 


V is of class ^provided a < 8(V) < m 


V is of class Cy? provided m < 8(F) < {3 . 


Further assume that V a and F b are vertices that are classified as c ^ and c p , respectively. 

If we now consider the values of 8 varying linearly along the edge joining V a and we 
could choose the mid-edge point to be the point where 8 becomes equal to m which is the 
point where the classification changes from c a to cp . That would be 


m 


m ab = 


8(V a ) 


W b +\ 


8(V b )-m 


.8{V b )~8{V a )) ^8(V b )-8(V a ), 


We have also used the following approach which is based upon the idea of a preference or 
probability matrix. The user specifies the off-diagonal values of an Mx M matrix 

^ = ( Py ) • These values serve as the weights for computing the mid-edge points. Let 
vertex V a be of class c a and V b be of class cp be vertices of the same tetrahedron. Then 
the edge joining V a and V b will intersect the separating surface at p ap V h + p Pa V a . Since 
the separating point must be a convex combination of the vertices we require that 
0 — Pij — 1 an d Pij + p jj = 1 , i ^ j . The interpretation of the matrix P can be in terms of 
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the "strength" of various classes relative to other classes, or it can be used to cause the 
separating surface to come close to (or stay further away from) certain classes of points. 
For example, it may be desirable to not overestimate the volume associated with a 
particular class, and in that case the values in the row of the matrix P associated with that 
class should be close to zero, forcing the separating surface close to the vertices of that 
class. 

2. Examples 

The first example has three regions. Points above the plane z = 0 , and outside the 
sphere x 2 + y 2 +z 2 =0.25, are of one type. Points below the plane and outside the 
sphere are in a second class and the points inside the sphere form the third class. Over the 
domain {( x,y,z ): - 0.025 < x < 0.625, - 0.625 < y < 0.625, - 0.625 < z < 0.625} we 
formed a grid of size 14 x 26 x 26 and classified the points on this grid according to the 
ideal. We then applied our algorithm, using the 5 tetrahedra per cube split. In this case 
we used the formulas of equation (1) for determining the mid-edge, mid-face and mid- 
tetrahedron points. The results are shown in Figure 7. 



Figure 7. An example with three regions. 
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One of the features of our algorithm is that it is designed for scattered data. Our next 
example illustrates its use in that context. Because the algorithm we used tetrahedrized 
the convex hull, points near the boundary may be in tetrahedra with large aspect ratios. 
This causes distortion around the boundaries, so in Figure 8 the separating triangles near 
the boundary have been deleted. We note, however, that if the proper tetrahedrization is 
performed, the separating triangles could be processed using subsets of the entire data set 
because our algorithm guarantees a proper match across tetrahedron boundaries. For 
Figure 8, we generated 2000 random points in the region 
{(x,y,z): - 0.2 < x < 1, - 1 < y < 1, - 1 < z < 1} . These points were then classified 
according to the same scheme as for the previous example. The point set was 
tetrahedrized and our algorithm applied with separating points being taken according to 
equation (1). To avoid the distraction of poor edge behavior, we then eliminated each 
triangle in the separating surface whose median point fell outside the region 

0 ^ x < 0.75, - 0.75 < y < 0.75, - 0.75 < z < 0.75} . The results are given in 
Figure 8. The separating surface is necessarily jagged, but the proper character is shown. 



Figure 8. An example with three regions, random points. 
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The final example has five different regions. These regions are defined relative to 
several conic surfaces, and the volumes are described sequentially, with a given class 

overriding a lower numbered one. Above the paraboloid z = 0.5(x 2 + y 2 ) the class is 1, 

while below (or on) the paraboloid the class is 3; inside the sphere 

x 2 + y 2 + (z - 0.75) 2 = 0.4 , the class is 2; inside the sphere x 2 + y 2 + (z + 1) 2 = 0.8 , the 

class is 4; and finally, inside the ellipsoid 2x 2 +(^-0.5) 2 +(z-0.l) 2 = 0.6 , the class is 
5. We formed a 41x41x41 grid over the domain 
{( x,y,z ): -l<y<l,-\<z<\} and classified the points according to the 

definitions of the various regions. Using the 6 tetrahedra per cube split, we ran our 
algorithm on this data using the P matrix P tj = l-0.2(y -/) for i<j. Because of the 
dense set of separating triangles, the results are shown as a shaded object in Figure 9. 



Figure 9. A surface separating five regions. 
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3. Summary and Remarks 

An algorithm which produces a solution to a problem similar to that discussed here is 
described in [1], This algorithm is based on case tables for the various configurations of 
classified vertices of a cube (or voxel) rather than a tetrahedron as used here. The number 
of equivalence classes of configurations (under rotations and possibly also mirroring 
transformations) for two and three classes of vertices on a cube is manageable but for more 
than three classes the authors mention that a table look-up approach is probably not viable 
due to the large number of different configurations. 

It is possible to use the techniques of this paper to build the case tables for cube 
scenario, but the results are different (possibly simpler) than what is presented in [1] due to 
the difference of using trilinear interpolation compared to piecewise linear interpolation 
over tetrahedron. For example, in the case where S = {( x,y,z): R.(x,y,z ) = R 0 (x,y,z)} is 
the separating surface with 

R.(x, y, z) = xyz + (1 - jc)(1 - y)z + (1 - x)y(l -z) + x(l - y)(l - z) and 

Ro(x,y,z) = (1 - x)(l - jv)(l - z) + xy(\ - z) + x(l - y)z + (1 - x)yz in the case of trilinear 
interpolation and R, and R 0 appropriately defined for the piecewise linear case, the 
piecewise linear case (using the 6 tetrahedra split) leads to a separating surface which is a 
manifold topologically equivalent to a plane, but the trilinear interpolation method leads to 
a surface that is topologically equivalent to the intersection of three planes! 

The algorithm presented here assumes that the data has been segmented into various 
classes and cannot be applied until this is accomplished. The problem of segmenting data 
is a highly nontrivial and currently unsolved problem. In no way does this present simple 
algorithm add to the solution of this problem, but possibly a more general algorithm which 
produces a tetrahedrized volume representation of the regions for different classes could 
be a useful tool in this regard. In a future paper, we will present such an algorithm. 
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